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Analysis of Field-Aided, Charge-Coupled 
Device Transfer 


By J. McKENNA and N. L. SCHRYER 
(Manuscript received August 7, 1974) 


We study the numerical solution of a nonlinear, partial-differential 
equation that describes charge transport in a model of a charge-coupled 
device (CCD). This model differs from previous models in that field-aiding 
of the transfer is taken into account. Although a derivation of the transport 
equation 1s given, the main emphasis in the paper is on the numerical 
techniques involved, and no actual numbers are presented. Actual numerical 
results based on the techniques developed here can be found in several 
recent design studies. The equation, which is parabolic, has one space 
dimension and one time dimension. Galerkin’s method, with standard 
chapeau functions, ts used to discretize in space. This results in a very 
stiff system of nonlinear, ordinary, differential equations. To solve these 
equations, we use a first-order backward Euler scheme coupled with 
extrapolation. A number of alternative schemes were tried and found to be 
inadequate. 


I]. INTRODUCTION 


In this paper, we study the numerical solution of a nonlinear, partial- 
differential equation that describes charge transport in a model of a 
charge-coupled device (ccp). The emphasis is on the numerical tech- 
niques involved, although a derivation of the equation is given. The 
reader is referred to other papers where the solutions are used in device 
theory and design.!? We briefly summarize the physical background 
of the equation first. 
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A knowledge of the dynamics of charge transfer in a cop is, of 
course, central to a complete understanding of its operation. A calcu- 
lation of the motion of charge in a ccp, starting from the coupled, 
nonlinear Poisson and charge-conservation equations and taking into 
account the full geometry of the device, has so far proved impossible. 
However, Strain and Schryer* and, independently, Kim and Lenz- 
linger* developed and studied an approximate, one-dimensional model 
of charge transfer in a ccp. The original analysis considered motion 
owing only to diffusion and the mutual repulsion of the charge carriers. 
lield-aided transfer was ignored. Since these original studies, a number 
of other authors have studied the effects of field-aiding.*-8 In Refs. 
5, 6, and 8, as in the original papers,?* an infinite sink for the charge 
at one end of a cell is assumed. The assumption of an infinite sink 
rules out charge “bunching,” which in certain situations is an im- 
portant effect (for an example of this, see Ref. 1, Fig. 8). In Ref. 7, 
the assumption of an infinite sink is not made. In this paper, we 
extend the original work** to include field-aiding and more realistic 
boundary conditions. Our model can describe both surface? ccps and 
buried-channel ccDs (BccDs). We do not include the effects of surface 
traps, since the main application! was to BccDs. We feel the numerical 
scheme described here has advantages over that used in Ref. 7, where 
essentially the same model as ours was used to study surface ccbDs, 
with the effect of traps included. Calculations using our methods 
show that Bccpds, which can be fabricated with present technology, 
should be extraordinarily fast and efficient and have reasonable 
charge-carrying capabilities. Transfer times of 1.8 ns are predicted for 
a two-phase device having 10-um-wide electrodes.! Slower but similar 
results are obtained for surface devices. 

Strain and Schryer* solved, by the method of finite differences, a 
transport equation quite similar to the one we study here. However, 
their method of solution proved inadequate when applied to our 
equation. It is possible to obtain solutions of the transport equation 
as follows. We use Galerkin’s method" with standard chapeau func- 
tions in space. We treat the time behavior by polynomial extrapolation 
to the limit of the results of a first-order, fully implicit (nonlinear), 
finite difference scheme. Although the equation only roughly models 
the true physical situation, an accurate knowledge of the solution as 
it varies over many orders of magnitude is necessary if it is to be of 
any use. This requirement makes the numerical solution of the equa- 
tion difficult. Many other schemes were tried, and the above method 
is the only one we found that could solve the problem. 

The equation of charge transport is derived in Section II, although 
some more complex details are given in Appendix A. The technique 
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for numerically solving the equation of charge motion is given in 
Section III, with some details in Appendix B. Questions of existence 
and accuracy are discussed in Section IV, along with the use of poly- 
nomial extrapolation. An outline of the theory of extrapolation is 
given in Appendix C. The method by which initial solutions are 
obtained is the subject of Section V. Finally, in Section VI we discuss 
several other schemes by which we tried to solve the equation of 
charge motion and which failed. 


Il. DERIVATION OF THE TRANSPORT EQUATION 


We refer the reader to the literature for a discussion of the principles 
of operation of either surface ccDs® or BccDs.'° Basically, however, both 
are devices that move packets of charge from under one electrode to 
under another electrode by suitably changing the voltage on the 
electrodes. 

As in Ref. 3, we assume that the charge can be described by a charge 
density q(x, t). Here, x is the distance under the plates (see Fig. 1) 
and ¢ is the time. Then, as we show in Appendix A, the component of 
the electric field along the direction of motion of the charge, which is 
due to the mutual repulsion of the charge, is 


Et a Sqz- (1) 


The elastance S is assumed to be a constant independent of x and ¢. 
In all that follows, we use subscripts to denote differentiation; thus, 
Qe = 0q(2, t)/dx, etc. Equation (1) holds for both surface and buried 
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Fig. 1—Schematic of a ccp showing relation to device of z-coordinate in transport 
equation. 
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channel devices, although the values of S are different in each case. 
Expressions for S are given in Appendix A in terms of the physical 
parameters of the devices. 

Let ¢(2, t) be the given driving potential due to the voltages applied 
to the electrodes. For a surface ccp, ¢ is the electric potential at the 
oxide semiconductor interface, while, for a BccpD, ¢ is the potential at 
the potential minimum of the buried channel. In most applications, 
we have approximated ¢ by the potential in the ccp in the absence of 
any mobile charge.??)8 

The total field along the direction of motion is 


EL, = — Sq: — @Qz. (2) 
The current density is? 


J (x, t) = ques = Daz, (3) 


where D is the diffusion constant and yu is the mobility, which we also 
assume to be constant. If we substitute (2) into (3) and make use of 
the Einstein relation D = (kT/e)y = ap, then 


J (x,t) = — ul(a + Sq)qz + Qee]. (4) 
If we substitute (4) into the charge-conservation equation,“ 
q: + Jz = 0, (5) 
we get the desired transport equation, 
qe = uL(a + Sq)ge + Q¢eJe- (6) 


The appropriate solution of (6) satisfies an arbitrarily given initial 
distribution of charge q(x, 0) and the boundary conditions J (0, #) 
= J(L,t) = 0. The boundary conditions state that there is no charge 
flow into or out of the device at either end. L is the length of the device. 

It is convenient to write (6) in terms of dimensionless quantities, 
as in Ref. 3. Let 


7 = t/(L?/pv), y= 2/L, w= S8q/vo, © = v/v, B=a/v%, (7) 
where vp is a reference voltage. Then (6) becomes 
w, = [(w + B)w, + wh,]y. (8) 


As it turns out, there seems to be no natural voltage unit in the 
problem (Ref. 3), so we typically pick vp = 1 volt. 

Physically, the quantity of interest is the total charge present 
between any two points 0 S y1 < y2 S 1. This suggests that, instead 
of w(y, 7), we consider 


Quy, 1) =f w(e, nae. (9) 
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If we integrate eq. (8) with respect to y from 0 to y and make use of 
the boundary condition J (0, ¢) = 0, we get 


Q, = (Qy ae B)Quy oF Qy Py. (10) 


Since the right-hand side of (10) is just proportional to J(y, 7), we 
see that Q,(1, 7) = 0. From this last remark and (9), it follows that 
the correct boundary conditions on Q(y, 7) are 


Q(0, 7) = 0, Q(1, 7) = Qr = const. (11) 


The appropriate initial condition is determined from w/(y,0) by 
setting 7 = 0 in (9). The transport problem we wish to solve is, 
thus, eq. (10), subject to boundary conditions (11) and given initial 
conditions. This is a much simpler problem than attempting to solve 
(8) for the charge density. 


Ill. SOLUTION OF THE TRANSPORT EQUATION 
We simplify the notation slightly by setting 


V(y, T) 7 b,(y, a); (12) 
and note that (10) can be written 


1 a F i 
— BQuy = 9 oy (Qy) = VQy oe Q, = 0. (13) 


If we multiply both sides of (13) by a continuous, piece-wise differ- 
entiable function f(y) which satisfies f(0) = f(1) = 0, integrate the 
result from 0 to 1, and integrate the terms containing second deriva- 
tives by parts, we obtain (letting f’ = df/dy) 


ih {L6Qy + 3(Q.)*If’(y) + [-¥Qy + Q-]f(y)}dy = 0. = (14) 


Equation (14) is the starting point for the application of Galerkin’s 
method, because any twice-differentiable function Q(y, 7) that 
satisfies (14) for all continuous, piece-wise differentiable f(y) satisfying 
f(0) = f(1) = O must also be a solution of (13). 

We now discretize in space by introducing a net {yi, ye, ---, yw} 
on [0,1] and a set of standard chapeau functions f;(y), 1 <j SN, 
as pictured in Fig. 2 and defined in Appendix B. In all that follows, the 





Fig. 2—Discretization of the space interval and the corresponding chapeau functions. 
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net {yi, ++, yv} is assumed to be given and fixed. In terms of the 
basic chapeau functions, we define approximations to the solution and 
external field: 


Oly, 2) = E Qilnfily) + Ont), (15) 
. N 
by, ) = Luni). (16) 


Note that Q(y, 7) has been constructed to satisfy the boundary 
conditions, O(0, r) = 0, O(1, r) = Qr. The functions Q;(r) are yet 
to be determined, but we require that they satisfy the initial conditions 


Q;(0) = Q(y;, 0). (17) 


Because of (17), O(y, 7) satisfies the correct initial conditions at the 
mesh points: O(y;,0) = Q(y;, 0). We define ¥;(7) = W(y;, 7), so that 
V(Yyi, T) = V(Yi, t). 7 

To determine the N — 2 functions Q;(7r), we require that QO(y, 7) 
satisfy (14) for each of the N — 2 choices of f(y), f(y) = fi(y), 
25 7<N —1, with y(y, 7) replaced by p(y, 7). This yields a system 
of N — 2 first-order, nonlinear, ordinary differential equations for 
the Q;(r). This technique has a robust history and has been applied, 
not only to many problems of the same type as (13), but to other 
types of problems as well. The idea is quite simple: Let the approxi- 
mate solution be a linear combination of the functions f;,2 <j S$ N—1 
and then make the left-hand side of (13) orthogonal to each of these 
functions. In geometrical terms, this means making the left-hand side 
of (13) orthogonal to the span of fo, ---, fy_1, denoted by (fe, ---, f—1), 
in £7[0, 1] in the usual inner product: (f, g) = Jo f(y)g(y)dy. Then, 


crudely speaking, as more points y; are chosen, (fo, ---, fy-1) spans 
more of £7(0, 1] and the left-hand side of (13) must go to zero as 
N —~, so long as it remains orthogonal to (f2, ---, fy—1). 


If we carry out the substitution of (15) and (16) into (14), with 
fy) =fity), 2 StS N — 1, the N — 2 equations result: 


306 f fifidy — Sve f° fiitedy + On f° fifetndy | 
N-1 1 i ae 1 
+45 iQ: ffilifay +X, Ose) fo fatedy 
= = |F fi sictavray + Oe (6 fl Sify 
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In (18) Q;(r) = dQ;/dr. The values of the integrals appearing in 
(18) are listed in Appendix B. If we define h; = yi — yii,2 SiS N, 
and substitute into (18) the values of the integrals, we get 


(hiQe—-1 + 2(hi + hig Qi + hiziQinr)/6 
+ Q:a[—B/hi + (bi-1 + 2p:)/6] + Q:[BC/hi + 1/hix1) 
+ (Wega — Wi-1)/6 + 63,n-107/ (hy)? ] + Qinal—B/hiss 
— (Wirr + 2:)/6] + 34 (Qi — Qi-1)?/h7 — (Qi4a — Q:)?/hi41} 
= 8;,v-11Q7/ (2h) + Qr(B/hw + Yw_-1/3 + Yw/6)], (19) 


where 6;,v-1 is the Kronecker delta function. These equations hold for 
2<isN —1if we let Qi(r) = Qn(r) = O. The nonlinear ordinary 
differential-equation initial-value problem given by (17) and (19) 
represents the spatial discretization of (13) and must now be solved 
for the Q;(7),2 Sj SN —-1. 
We use a fully implicit finite difference scheme in time (backward 
Euler). Let 
Q? = Q;(nAr) (20) 
for some choice of Ar > 0. We then let Q;(nA7r) be approximated by 
(Q7*? — Q?)/Ar and set Q; = Q?t!, y; = v7t in (19). On rearranging, 
we obtain the fully implicit, first-order, finite-difference scheme for 
solving (19) in time: 
THQ + THQI + TROT + AQ?! — ory 
— Ai (QHi — Q7*1)? = RIM, (21) 


where 
Tit! = — B/hy + (Witt + 2p27)/6 + hi/ (647), (22a) 
Tit) = B(1/hs + 1/higs) + WHA — WEET)/6 + (hi + higt)/ (847) 
+ 6:,n-1Q7/(hy)*?, (22b) 
Tie == B/hiss — (watt 7h Q2y7tt) /6 fe hi+1/(6Ar), (22c) 
A; = 1/(2h), (22d) 
R?t! = 6; y_i1[Q3/(2hx) + Qr(B/hy + Prt /3 + vrt/6)] 
+ (hiQ?-1 + 2(hi + hig. Q? + higrQ?-1)/(6Ar).  (22e) 
Equations (21) hold for 2 $i N —1, n=0, 1, 2, ---, with the 
assumption that Qj = Qy = 0, n = 0, 1, 2, --- and with the initial 


conditions Q3 = Q(y:i,0),2 Sis N-—1. 

We now find the solution of the nonlinear system of eqs. (21) for 
fixed n by an iterative Newton method. We drop the superscript n 
denoting the time step, and for fixed n denote by Q:(m),2 S75 N —1, 
the mth iterate of the solution of (21). To obtain Q;(m + 1) from Q;(m), 
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we set Q;(m + 1) = Q:(m) + ri(m), substitute this into (21), and 
linearize the resulting equations for the 7;(m) : 


{Tir — 2A, [Qi(m) — Qi-1(m) ]}7s-1(m) 
+ {Tis — 2A igs [Qir1(m) — Q:(m) ]}rin1(m) 
+{Ti2 + 2A.:(Qi(m) — Qi-1(m) ] 
+ 2A i4i[Qizi(m) — Qi(m) ]}ri(m) 
= Ri — {TaQ-a(m) + TieQi(m) + Ti3Qi41(m) 
+ ALQi(m) — Qer(m) P — AtalQir(m) — Qi(m) PF}. (28) 


These equations hold for 2 $7 S$ N — 1 with ri = ry = O. This is a 
tridiagonal system of linear equations. Reference 15 contains a concise 
analysis and very efficient method of solution for such a system of 
tridiagonal equations. 

In practice, the initial estimate of the solution Q?t! to (21) is taken 
to be Q? from the previous time step. So, if Az is chosen sufficiently 
small, the Newton sequence generated by (23) should converge and 
do so quickly. 

What we have described so far is a method for discretizing (9) and 
(10) in space and time, giving the nonlinear system of eqs. (21), and 
we have proposed an iterative scheme, given in (23), for solving (21) 
at each time step. In the next section, we study the feasibility and 
accuracy of the method. 


IV. EXISTENCE AND ACCURACY 


We shall show that iteration (23) can be carried out as long as the 
following conditions are satisfied: 


O0OSQ508--- oe mS Oil 2 ty, 1124) 


AAC Ol S5%. 2SvsN, (25) 
Cyi-1,.yi ny FO, « 

These conditions are sufficient to ensure the existence of a solution of 
eqs. (23) for each n. We have not proved it, but in practice they also 
seem to be necessary. These conditions do not show that the iteration 
(23) must converge, merely that it is well defined. In fact, if the initial 
estimate of the solution of (21) is too far off, then in practice the 
Newton sequence given by (23) may well not converge, and it is 
necessary to choose Ar smaller so that Q? provides a better estimate of 
OFe. 

The monotonicity condition (24) on Q? is merely a necessary conse- 
quence of the definition (9) of Q?, since w(é, 7) 2 0 by definition. The 
mesh restriction (25), however, is apparently new and fundamental. 
In practice, if (25) is violated, even at only one point and by a ‘“‘small’’ 
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amount, the solution produced, if any, is highly erratic and non- 
monotone, and may even be negative. 

We now prove that conditions (24) and (25) imply that the matrix 
of eqs. (23) is strictly diagonally dominant.!* From this, we can con- 
clude that the matrix has an inverse,!* so the equations have a solution. 
From (22a) to (22c) and (25), we see that 


Tia t+ Tie + Tis = (hi + hitt)/(247) > 0, (26) 
and 
Tin 2 (hi + higi)/(847) > 0; 
La = h;/(6Ar7), Tis S hisi1/(6Ar). 


(27) 
Because of the monotonicity property (24) and the fact that Ts. > 0, 
it is easy to show that 

AT; = Ti —|Tal|—|Tis| >0 (28) 
implies the diagonal dominance of (28): 


[Tis + 2A,[Q:(m) — Qi-1(m)] + 2A 1[Qi41(m) — Q:(m) ]| 
> | T1—2A,[Q;(m) = Q:-1(m) ]| 
+ |Tis — 2AsilQi+1(m) — Qi(m)]|. (29) 
To show that (28) is true, we consider the four possible sign combina- 
tions of 7; and 73 and use (26) and (27): 
(t) Ti: > 0, Tis > O. 


AT; = Ti — Tia — Tis = (he + higit)/(247) 

— 2(7T 1 + Tis) 2 (hi + higi)/(6Ar) > 0. 
(7) Tia > 0, Tis < 0. 
AT; = Tig — Ta + Tis = (h; + hizs)/ (247) 


hi hiss 
= oe Se ee 
22a 2 Gh Ap 





(412) Tia <0, Tis > 0. 
AT; = Tig + Tar — Tis = (hi + hits) / (247) 


es es 
20s 2 2Ar 


Hie Meets 


6Ar ; 








(2v) Tin < 0, T 53 < 0. 
AT; = Tig + Ta + Tis = (hi + Aigi)/(247) > 0. 


This completes the proof of the diagonal dominance of (28). 

We now discuss the accuracy of the spatial and time discretizations. 
It is well known (see Ref. 11) that the Galerkin procedure, using 
chapeau functions, is accurate to 0(h?), where h = max,h, and O(h?)/h? 
represents roughly an upper bound on Q,,(y, 7) over [0, 1] X [0, ~). 
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We shall not go into the proof of such results here. Rather, a heuristic 
but useful analysis of the error is presented. The 0(h?) accuracy, ba- 
sically, comes from the fact that replacing Q(y, 7) by its interpolant, 


© Qs fy) + Onfv(y), 


results in such a O(h?) error by using Taylor’s theorem on each of the 
intervals [y:, y:41], 7 = 1, ---, MN —1. A similar statement can be 
made about y(y, 7) and its interpolant. For the sake of clarity, assume 
that the mesh is uniform with h; = h, i = 2, ---, N. Then standard 
finite difference arguments show that (18) is a spatial finite difference 
approximation to a function Q*(y, 7) obeying 


Q; = Q(B + Q) — ¥Q, + 0(), (30) 


where 0 involves terms of the form Q}, and its higher-order derivatives, 
a”*"/ay"d7". Then, intuitively speaking, since Q(y, 7) solves (30) to 
within 0(h?) and Q*(y, 0) — Q(y, 0) = O(h?), we must have Q*(y, 7) 
— Q(y, 7) = OCA’). 

Even though (30) is based on the assumption that the spatial mesh 
is uniform, it shows clearly that the h; must be small in any region 
where any of the derivatives (0"+"/dy"dr")Q,, are large. Physically, 
such regions are precisely those regions where the field y(y, 7) is large. 
This makes restriction (25) quite reasonable, since (25) requires a 
smaller spatial mesh where the field y is large. In fact, we can estimate 
the number of points V,, required by (25), using a variable mesh, in a 
potential rise of v volts: (25) requires that y change by no more than 
28 = 1/20 (at room temperature) over any mesh interval. Then, for 
example, a potential rise of 5 volts will have 100 points y; modeling 
it. So (25) itself forces a fairly accurate representation of y and hence, 
indirectly, of Q. 

However, the time mesh is another matter altogether. The time 
difference scheme is only first-order accurate and the local time 
behavior of Q near large values of y is rather bad. Thus, application 
of (21) to (23) alone to solve the problem gives rather poor results. 
For this reason, we have used polynomial extrapolation to the limit of 
the results of the first-order scheme (23). A brief discussion of the 
extrapolation process is given in Appendix C. Ironically, polynomial 
extrapolation was used because rational extrapolation converged so 
quickly to the solution that it led to very large Az choices (see Ref. 17 
for the Ar monitoring mechanism) which, in turn, led to iteration (23) 
not converging or taking a very long time doing it. So, even though 
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polynomial extrapolation is “slower” than rational, it is “better” for 
our purpose here. 


V. CALCULATION OF Q(y, ~) 


In most cases of interest, the initial condition for (10) is chosen 
as an equilibrium solution Q(y, ©) corresponding to a time-inde- 
pendent potential @(y). It is convenient in these cases to solve for the 
corresponding w(y, ©) = w(y) and then integrate to get Q(y, ©). 

Setting w, = 0 in (8) yields 0 = [(w + 8)w, + w@, ],, which, when 
integrated twice from 0 to y with the aid of the boundary condition 
J(0, ©) = 0, yields 


F(w) =w+Blnge + ¥(y) = 0 (31) 


for some constant Co. Let yo be any point in [0, 1] such that w(yo) > 0. 
Then 


Oy ws) exp ( Ba) Feta) ): (32) 


Thus, given ®(y) and a single value of w(yo) > 0, the entire equi- 
librium distribution w(y) is determined. Note that w(y) > 0 whenever 
@(y) is finite. 

To find w(y) from (31) we use Newton’s method. An initial guess at 
the solution w(y) > 0 is made. The solution is then iterated, the 
(n + 1)th iterate being related to the nth by 


wey) | + Seg = 3) +8 [I In (“oP i 


Since F’(w) = 1+ 8B/w > 0 and F’’(w) = — B/w? < 0, we see that 
F(w) is a concave, monotone-increasing function for w > 0. Thus, the 
Newton sequence generated by (33) will converge to the solution (31) 
no matter what initial w(y) > 0 is chosen. 

Once the w(y;), y: in the Galerkin net {yi, ---, yv} are found using 
(33), Q(y:) may be found by the trapezoidal rule for integration. This 
is consistent with the representation of Q by the chapeau functions, Q, 
since the trapezoidal rule is exact for chapeau functions. 


VI. ATTEMPTS THAT FAILED 


The first attempt at solving (10) was via the finite difference scheme 
of Ref. 3. It was impractical because the spatial mesh restriction (25) 
appeared there, also, forcing the spatial mesh to be very small in some 
regions, although it could be quite large in others. Since any non- 
uniformity of mesh size in a central finite difference scheme leads to 
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only first-order accuracy, we were then left with a very fine mesh over 
the entire interval [0, 1]. This required tens of thousands of points in 
the spatial mesh, far too many to be practical. 

After going to Galerkin’s method in space, which has second-order 
accuracy even with a nonuniform mesh, the solution of (19) posed 
another problem: It is an extremely ‘‘stiff’’? system of ordinary differ- 
ential equations, with the coefficients A; ranging typically from 10! 
to 10", This is a reflection of the locally quick time and spatial changes 
in Q(y, 7) when y is large, this fact being transmitted to the h; by (25). 
For this reason, any attempt to linearize (19) between time steps for a 
finite difference scheme in time led to failure—the solution is nowhere 
near linear over reasonable time intervals when y is large. The symptom 
of this problem, in practice, was that the Az required in the polynomial 
or rational extrapolation process for these linearized schemes was 
extraordinarily small, requiring in one case more than 10" time steps 
to cope with a single 5-volt potential swing. 

Once a nonlinear approach to the solution of (19) was recognized 
as probably the only route left, the most obvious ‘‘accurate’’ scheme 
to use is a fully nonlinear Crank-Nicholson solution of (19). A small 
digression on this scheme in a simple case is useful here. For the linear 
system of ordinary differential equations, 


u’ = Au, (34) 


where u is a vector and A a matrix, the Crank-Nicholson approximation 
to the true solution, u = e47Us, is 


u(nAr) & (I + $AAr)"(I — $AAr)~"Uo. 
This is based on the approximation™® 
eA4t = (I + $AAzr)(I — 3AAz)71. (35) 


Letting u(nAr) = (ui, ---, ux)”, this corresponds to the standard 
finite difference formulation of the Crank-Nicholson scheme: 


(ug? — uf)/Ar = 3(Au™4 + Au), 17 SN. 


A nonlinear generalization of the above scheme for (19) would have 
an error of the form C(Ar)?; however, C is very large. This is most 
easily seen by considering (35) for real AAz very large (positive or 
negative). That relation then states that e44” = — 1, which is an 
exceedingly bad approximation. For a “stiff” system, (34) [or (19) ], 
one that has a wide spread in its eigenvalues for A, the above reasoning 
indicates that the Crank-Nicholson scheme would give very poor 
results unless Ar is very small. In practice, as before, the symptom of 
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this problem was very small Ar choices by the extrapolation routines— 
the same problem that would have required 10" time steps in a linear 
scheme would have required ‘‘only” 108 with Crank-Nicholson. (In 
this matter, see also Ref. 19.) 

In all, more than 12 different schemes were programmed and tested 
on this problem, (9) and (10), with the result that only the one de- 
scribed in Sections II to V is effective for the wide range of y distri- 
butions required to model both surface and buried-channel ccDs. 
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APPENDIX A 


In this appendix, we derive eq. (1), the fundamental equation of 
the Strain-Schryer model, for the case of a Bccp. We choose coordinates 
as shown in Fig. 3; the z-axis is parallel to the oxide-semiconductor 
interface and the z-axis is directed into the semiconductor. The 
potential in the oxide is ¢o(x, z) and the potential in the semiconductor 
is gi(a,z). The permittivity of the oxide is e.:, that of the semi- 
conductor e,, and the thickness of the oxide is 6. 

In the special case where all the properties of the Bccp are inde- 
pendent of x, the potential in the presence of the inserted charge q 
has been calculated by Kent?" and Schryer.”4 They showed that the 
value of the potential at its minimum in the buried channel is approxi- 
mately a linear function of the charge q, ¢1 = Sog + Vo, for all values 





Zz 


Fig. 3—Coordinate system involved in calculating the potential of a line charge. 
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of q in the operating range of the device. The elastance Sp and V» are 
independent of q but depend on the oxide thickness 6 and the semi- 
conductor doping, and V> also depends strongly on the electrode 
potential. 

In the general case, the Strain-Schryer model assumes that the 
field in the x-direction in the channel can be approximated by the sum 
of two terms. The first term is obtained from the above expression for 
gi by assuming q a function of x and differentiating, 


o 
Ex = — So = (36) 


The second term takes into account the field at x resulting from the 
charge at other points x’ in the channel. Because of the metallic elec- 
trodes, the charge at x’ will induce image charges that will tend to 
shield the field at x. For this purpose, we first calculate the potential 
of a unit line charge located at « = 0, 2 = 7 > 6 in the semiconductor. 
The plane z = 0 is assumed to be a perfect conductor at zero potential, 
and the oxide and semiconductor are assumed uniform. We can write 
down a solution of Laplace’s equation in the form 
cole, 232) = [ r(a) soa eed, O52 58, (37) 


gi(X, 259) = — i W(x, 237) + ie S(a)erlal 2-9 taxa 
6S2< 0, (88) 
where 
W(x, 230) = In {a2 + (2 — 0)9} — In {a2 + (2 +0)*}. (39) 


The function V(a, 2; 7) has the correct singular behavior at x = 0, 
z = n and is harmonic everywhere else in -—-~ <x4< w,652< ~, 
The boundary condition g(x, 0; 7) = 0 is satisfied, and the unknown 
functions r(a) and s(a) must be chosen so that the boundary conditions 


rs) 0 
go(%, 330) = gilt, 331), €or a (#, 830) = & 3 (2,850) (40) 
Oz Oz 
are satisfied. It is straightforward to show that 


W(x,6;y) = — 2 e—nlel soualt eda, (41) 


ae | 
‘- (z,5;9) = —2 / e-"l@l cosh |a| eda. (42) 
If we substitute (37), (38), (41), and (42) into (40), and Fourier- 
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transform with respect to x, we obtain two linear equations for r(a) 
and s(a). The solution of these two equations yields 


r(a) = te Gr-Dlel} (eqs te )elel® 4+ (eae — «Je !al2}—1, (43) 


sinh|a|6 | 


Te 





(44) 


On substituting (43) and (44) into (87) and (88), we obtain the 
desired result. If we expand r(a) and s(a) in powers of e~!#!®, the 
Fourier integrals can be evaluated, and we can express the potential 
as the potential resulting from an infinite array of image charges. Since 
this result is not needed, we do not give it here. 

In the buried-channel case, we need the potential resulting from a 
two-dimensional charge distribution. Let the density of this distribu- 
tion be p(é 7). Then g(x) = S p(x, n)dy is the charge appearing in 
eq. (36). Since the potential resulting from the image charges induced 
by a line charge at (& 7) in the semiconductor is g(a — & 237), we 
can now write down the second term of the field in the channel as 





Ba=— f [2 @- & zine ndter (45) 
From (38) and (89), 
— (@ — & 23) 
2S ait | e—& _ wot | 
2m | = Ee en)? eae 8)? af 2 eg) 


44 i as(a)enlal (e-Beiae-Hdy, (46) 


Since (0¢1/0x)(x — £,2;n) is singular at = 2, » =z, the main 
contribution to the integral in (45) occurs at this point. We expand 
p(n) in a Taylor series about x, keep only the linear terms in the 
expansion, and extend the limits of the & integral from — to . 
Since (0¢1/0x)(x — £237) is an odd function of x — &, the term 
involving p(x, 7) vanishes. A straightforward calculation shows that 
the ioe term is 


1 1) 99, 
~ ck feta le- alee adr-5(2-2)5t an 
The first integral can be transformed by the mean value theorem: 
S (@+— |e — a])e(2, adn = (2 +4 — |2 — alga), where 7 is a 


point in the interval of integration. In many cases, it is reasonable to 
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replace the factor z + 4 — |z — 4| by a constant 21, independent of z. 
For such cases, we have 


Pps {if +0(2-2)] 99. (48) 


€or Es Ox 
If we combine (36) and (48) we obtain (1), where 
S = Sot (1 — 6)/es + 6/ Ecc. (49) 


Here So must be obtained from a one-dimensional charge-insertion 
calculation,?*:4 and / must be estimated from the above formulas. 

It should be noted that, if we let p(é, 7) = p(&)D(n — 4) in the 
previous derivation, where D(x) is the Dirac delta function, and set 
y = 6, we should get the result of Ref. 3 for a surface device. However, 
in this case, (47) yields 6/e. for the correction term, while in Ref. 3 
the correction term is 26/(€, + €oz) [eq. (4) ]. This is because, in Ref. 3, 
in the expansion of the field in terms of image charges, only the first 
image was taken into account. 


APPENDIX B 


In this appendix, we list several results concerning the chapeau 
functions f;(y) : 


fy =0, OSyS yA, 
= (y — yj-1)/hy;, Yj 
= (Yin — Y)/Aj, Ys Yi4t, 
= 0, Yo41 = y = 1, (50) 


IA IIA 
< 
IA IIA 


where h; = y; — yj-1.- 
We list here a number of elementary integrals that are needed in 
obtaining eqs. (19) from eqs. (18). 


fo Gitdy = ty + Wher (51) 
[fifty = — Ahn (52) 
[fo Gp'dy = Chins + 19/3, (53) 
[ ddissdy = hiss/6, (54) 
fo Gdidy = ay* = nad (56) 
[Dina = Gan, (56) 
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[ G%iady = — 0), (57) 


[ G*fay = 0, (58) 
[ Gu¥¥idy = - 3, (59) 
[ ddidaady = 8, (60) 
[ DYndy = 3, (61) 
[ Sihiastinady = 3 (62) 


In all these expressions, f; = df;/dy. 


APPENDIX C 


In this appendix, we give a brief description of the extrapolation 
method for solving eqs. (19) in time. We used a linearized, backward 
Euler method for solving (19) in time. It is first-order accurate. That 
is, by using a time step of At to go from éo to t1 = ty) + mAt, the result- 
ing error at t; is O(At). See either Ref. 22 or Ref. 23 for the proof of 
such results. 

However, much much more is known about these methods. In fact, 
Stetter™* has shown, in a very general setting, that processes such as 
the above backward Euler technique give rise to expansions of the 
form 


T(At) = TCO) + ¥ 2)(Ad), (63) 
7=1 
where, for our problem, T(A?) is the value of the vector (QT, ---, QN)’, 


which is the value of our approximate solution at ti = t) + mAt, and 
the +; are vectors that depend only upon ¢ and 4. Thus, as 
At = (ti — to)/m goes to zero or, equivalently, as m goes to infinity, 
T(Az) not only converges, with error 0(At), to the true solution at 1, 
namely, T(0), but each component of T(A#) looks more and more like a 
polynomial in At. The process of extrapolation consists of simply 
computing several values, T(At), T(At/2), ---, T(At/p), and then 
passing a polynomial of degree » — 1 through these data points 
corresponding to each component. The value of these interpolating 
polynomials at the origin is the solution T(0), plus terms of order 
(At)”. Here p= is called the level of extrapolation. 
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By using polynomial extrapolation to the limit of the result of the 
first-order scheme (21), we generate a process that has an error of 
O[(At)?] when p levels of extrapolation are used. This extrapolation 
process is very well described in Ref. 25, and its application to the 
numerical solution of ordinary differential equations is also very well 
described in Ref. 17. It must be stressed that the underlying process, 
Gragg’s modified midpoint rule, which Bulirsch and Stoer extrapolate 
in Ref. 17, is not the one we are proposing to extrapolate here. That 
rule is second-order accurate and is actually unstable if the equations 
being solved are stiff. The first-order, linearized, backward Euler 
method we use here is highly stable under extrapolation, even for very 
stiff systems like (19). So Ref. 17 should be read with an eye to using 
extrapolation in solving ordinary differential equations and not to 
those peculiarities that Bulirsch and Stoer introduce to take special 
advantage of the nice properties of Gragg’s modified midpoint rule. 
The same technique we have used here to solve (13) was used in Ref. 26 
to solve a similar system. It is of interest that, for both these problems, 
polynomial extrapolation was found to be 15 to 20 percent faster than 
rational extrapolation. This is in contrast to the finding in Ref. 17 
that rational extrapolating is usually the faster of the two, at least 
when extrapolating Gragg’s modified midpoint rule. 
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Aluminum Oxide/Silicon Dioxide, Double- 
Insulator, MOS Structure 
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A double-insulator structure consisting of 500 A of vapor-deposited 
Al.O3 and 1000 A of thermally grown SiOz is used as the gate dielectric 
in a beam-lead-compatible, p-channel, MOSFET, silicon-integrated-circuitt 
technology. The AlzO; layer, 1n addition to serving as a sodium barrier 
and thereby providing a self-passiwated technology, results in a positive 
flatband voltage shift when compared to an S102 structure. The mechanism 
for this flatband voltage shift is the subject of this paper. 

The major experimental results obtained are (t) a negative charge 
exists near the Al,O3/StO2 interface, its magnitude being independent of 
the Al,O3 thickness but inversely proportional to the SiO, thickness, 
(it) the magnitude of the SiO2./St interface charge is inversely propor- 
tional to the SiO thickness, and (itt) a potential jump of about 1.26 volts in 
flatband voltage 1s associated with the addition of the Al2O3 layer. 

A physical model 1s proposed which assumes the existence of a constant 
voltage drop across the S10, layer during the Al,O; deposition and a 
corresponding charge buildup at the St02/Al,03 interface. 


1. INTRODUCTION 


The threshold voltage of an insulated-gate, field-effect transistor is 
directly dependent upon the properties of the gate insulator. A double- 
dielectric gate structure consisting of nominally 500 A of vapor- 
deposited Al.O; and 1000 A of thermally grown SiOz is the basis of a 
beam-lead-compatible, p-channel, MOSFET, silicon-integrated-circuit 
technology.1-* The Al.O; layer serves two functions. First, it is a 
diffusion barrier for light ions, such as sodium, and thus provides a 
self-passivated technology. Second, the Al,O; layer shifts the threshold 
voltage of the mos transistor in the positive direction (due to a flat- 
band voltage shift). For example, for a (100) oriented, n-type, 10- 
ohm-cm, silicon substrate, a flatband voltage of 0.0 volt is obtained 
with the dual-dielectric structure and a titanium metal gate, charac- 
teristic of the beam-lead metallization system, whereas with just an 
SiO. structure the flatband voltage is —0.8 volt. The more positive 
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flatband voltage capability provided by the Al2,O3 layer implies that 
MOSFET integrated circuits can be fabricated which have low power- 
dissipation properties and which are more easily interfaced with bipolar 
circuits. 

Many techniques have been used for the deposition of Al.O; films 
intended for application in an integrated-circuit technology. All of 
our considerations are restricted to Al,O3 films deposited at 900°C 
from an AIC]; source, the technique reported by Tung and Caffrey.! 
A brief review is given in Ref. 4 of the other Al.O3 deposition techniques 
that have been reported and the electrical characteristics of the films 
obtained. 

The electrical properties of the Al2O3/SiOe, dual-dielectric, gate 
insulator are quantitatively described in this paper; and, in particular, 
the mechanisms are delineated which cause the positive shift in flat- 
band voltage. The experimental approach was to do a parametric 
study of the flatband voltage of the dual-dielectric Mos structure, 
the two parameters of interest being the thicknesses of the Al.O3 
and SiOz layers. 

Two major conclusions were obtained from the parametric study. 
First, a net negative charge exists near the Al,03/SiO2 interface, and 
the magnitude of the charge is independent of Al,O; thickness over the 
range studied, but inversely proportional to SiO» thickness. Second, 
the magnitude of the normal interface charge associated with the 
Si0./Si interface has a component that is inversely proportional to 
SiO. thickness. 

A model explaining the origin of the negative charge at the Al.O3/ 
Si02 interface was developed, based on the assumption that the elec- 
trical conductivity of Al,O; at high temperatures (>3800°C) is much 
greater than that of SiO». As a result, during the high-temperature 
deposition of Al.O3, an electric field exists in the SiO. due to the 
Si/Al,03 contact potential difference, and the negative charge at the 
Al,03/SiO»s interface terminates this field. 

Recently, Aboaf, Kerr, and Bassous also reported the existence of a 
negative charge at the Al,O3/SiO: interface with the magnitude of the 
charge being independent of the Al.O; thickness and inversely propor- 
tional to the SiO, thickness. This is consistent with the insulator 
interface charge origin model we proposed’ and implies that this model 
has general applicability in dual-dielectric structures, since they used 
three Al,O; deposition techniques, all different from the technique 
used to obtain the Al.O; films studied in this paper. 

The organization of the paper is as follows. A simple flatband theory 
for dielectric structures is presented in Section II and space-charge 
formation in insulators is discussed in Section III. A theoretical 


688 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1975 


discussion of the Al,O3/SiO2 structure is given in Section IV and the 
experimental results are presented in Section V. A summary is given 
in Section VI. 


il. FLATBAND CALCULATIONS 


The flatband voltage of an Mos structure is defined as that voltage 
which must be applied to the metal electrode to produce a zero space 
charge or flatband condition in the semiconductor, and it is determined 
by the net charge density existing in the insulator system and the 
various interfacial barrier energies. To calculate the flatband voltage, 
the voltage across the insulators under flatband conditions is calcu- 
lated from the net charge density using Gauss’ law, and to this is 
added the voltage contributed by the various barrier energies. 

Consider the band diagrams of the metal/SiO2/Si and metal/Al.03/ 
8102/81 systems shown in Figs. 1 and 2, respectively. The various 
barrier energies for these systems are defined in the figures. We shall 
assume that at the SiO./Si interface in both structures there is an 
interface charge layer Q,,. The net charge density in the bulk of the 
SiOz. is assumed to be zero and the charge density in the Al,Q3 is 
denoted by pa (x). Using S.I. to denote the single-insulator system and 
D.I. the double-insulator system, it follows that the voltages across 
the insulators, V;, due to the charge densities can be written as: 





V,(S8.1.) = (Qee/ €or) T oz (1) 
ie! TAtT ox d. x 
V; (D. ve )= Qi |= ao | sia / me i pa(x’)dax’ 
Toxt+TaA qT a5 
= = Qss | + Eo + | ae i. Fat 4a? 54 (x)ade, (2) 


where (€oz, €4), (T'oz, Ta) = the dielectric constants and thicknesses 
of the SiO2 and Al.Os, respectively. In particular, €o2 = 3.9; «4 = 9.0. 

The applied voltage difference between the metal and the semi- 
conductor, the flatband voltage Vrg, for both structures can be 
written as: 


Vra(8.L.) => (Qss/ €or) T' ox + (@mioz —_ dB — by) (3) 
Vera(D.1.) = — VDL) + (¢n,4 + bi — GB — $y). (4) 


The insulator-insulator barrier ¢;; is assumed to be positive if it is as 
shown in Fig. 2. 

Consider the contact potential terms in eqs. (3) and (4). If we let 
W and X with the appropriate subscripts denote the vacuum work 
functions and electron affinities, respectively, of the various materials, 
and, if we assume that the energy band matching between two different 


DOUBLE-INSULATOR MOS_ 689 
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——Si0, 


——Si 


(a) 


SILICON SiOg METAL 








SILICON Si02 METAL 
(c) 
Fig. 1—(a) Cross-sectional view of the metal/SiO2/Si capacitor structure. (b) 


Band diagram associated with the structure. (c) Plot of assumed charge density 
present in the structure. 


materials is determined entirely by the difference in vacuum work 
functions, that is, 
Pm, ox = Wm == Kgs and op = X, — Xozy (5a) 
then 
dm,oz — OB — OF = (Wm ~ Xox) ~ (Xs _ Xox) — oF 
= Wm — Xs — $7, (5b) 
dm,A + bi — OB — Of = (Wm — Xa) + (Xa — Xoz) — (Xs — Xox) — Hf 
= Wm — Xs — Of; (6) 
or 
dmor — dB — bf = Gm,a + biti — GB — b7 = Wm — X, — Gy. (7) 


Note that under this assumption, the contact potential terms are 
independent of the electron affinities of the insulators and dependent 
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SILICON Si02 Aln03 METAL 





0 Tox Tox + Ta 
SILICON SiO, Aly 03 METAL 


(c) 


Fig. 2—(a) Cross-sectional view of the metal/Al.0;/Si0O2/Si capacitor structure. 
(b) Band diagram associated with the structure. (c) Plot of assumed charge density 
present in the structure. 


only on the difference between the metal and semiconductor work 
functions. Thus, if the work-function assumption is valid, an mos 
system involving a given metal and a semiconductor will have a 
constant flatband voltage after correction for Q;, independent of the 
number or nature of the insulators, unless space charge exists in the 
insulators. In other words, if the assumption is valid, a measured 
Vrg that changes when the insulators are changed implies space 
charge exists in the insulators. 


Ill. SPACE-CHARGE FORMATION IN INSULATORS 


In the previous section, the flatband voltage of a dual-dielectric 
Mos structure was calculated assuming a given distribution of space 
charge. The purpose of this section is to discuss one possible model 
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for the origin and spatial location of space charge in insulators, namely, 
the one we feel represents the most likely explanation for much of the 
space charge in the $i02/A1.0; system. In this discussion, we first 
summarize the proposal for space-charge formation suggested by 
Simmons for a metal-insulator-metal (mm) system.® Then this model 
is extended and applied to double-insulator systems, in particular those 
using SiOz and Al.O3. 

Suppose we form the mim system shown in Fig. 3a at sufficiently 
low temperatures that no charge transport occurs within the insulator 
and no charge exchange occurs between the insulator and the metallic 
contacts in a time period comparable to the experimental observation 
time. In this case, no space charge can form in the insulator because 
thermal equilibrium will not exist and the potential versus position 
will look as shown in Fig. 3b, where the insulator is represented essen- 
tially as a wideband insulator with conduction and valence band 
edges. Now assume that the system is heated to a sufficiently high 
temperature that charge transport can occur within the insulator and 
charge exchange can occur between the insulator and the metallic 
contacts in a time period that is short compared to experimental 
observation. In this case, thermal equilibrium will be established and 
there will be two extreme possibilities between which the system will 
equilibrate: either the characteristic length corresponding to a space- 
charge region at thermal equilibrium in the insulator, the electrostatic 
screening, or Debye length will be large compared to the insulator 
thickness, and the potential versus position will be virtually identical 
to that shown in Fig. 3b; or the Debye length in the insulator will be 
small compared to the insulator thickness, and such space-charge 
regions as shown in Fig. 3c will form near the two metal interfaces. 
In the latter alternative, a well-defined ‘J*ermi level’ will exist in the 
bulk of the insulator, as shown in Figure 3c, which will coincide in 
energy with the Fermi level in the metallic contacts in much the same 
way that Fermi levels coincide in a conventional Schottky barrier on a 
semiconductor. Clearly an Mim system in which the Debye length is 
short compared to the insulator thickness at any given time will 
always lie somewhere between the extremes indicated in Figs. 3b and 
3c, depending on thermal history, so that in such an insulator, space- 
charge regions will always exist in the vicinity of the metallic contacts. 
The magnitude of the charge will depend on the difference in the work 
function of the metal and the insulator and on the degree of thermal 
equilibrium which has been established. 

In the above discussion, it has been assumed that there is no net 
voltage difference across the mim structure. Very similar arguments 
can be presented for the case where a finite voltage exists between 
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Fig. 3—(a) Cross-sectional view of a metal/insulator/metal capacitor structure. 
(b) Band diagram of the sturcture if the insulator Debye length is assumed to be 
much greater than the thickness of the insulator. (c) Band diagram of the structure 
if ae insulator Debye length is assumed to be smaller than the thickness of the 
insulator. 


the two metal contacts; the voltage is either applied or is due to 
differences in the two involved metal/insulator potential barrier 
heights. Assuming an insulator which forms space-charge regions 
that are narrow compared to the insulator thickness, the initial 
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potential diagram is shown dotted in Fig. 4a and the final steady-state 
situation is shown as the solid lines. The only electric field required in 
the bulk of the insulator under steady-state conditions is the ohmic 
field associated with any current injected from the electrodes, typically 
a negligible field. 

Now suppose that the voltage is reduced as shown at ¢ = 0 in Fig. 
4b. Initially, the same space-charge that exists under applied voltage 
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Fig. 4—(a) Band diagram of a metal/insulator/metal structure, whose insulator 
Debye length is less than the insulator thickness, depicting the immediate and 
equilibrium band structure in response to the application of an external voltage. (b) 
Band diagram of the same structure depicting the time response of the bands after 
the applied voltage is reduced to zero. 
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will remain. After a time ¢; at elevated temperature, there will be a 
redistribution of charge as shown, the total amount of charge remaining 
fixed. Finally, if the temperature is further raised or after an additional 
time é2, carriers will be injected from one or both electrodes to bring 
the system to the equilibrium of Fig. 3c. 

The above argument for space-charge formation in insulators is an 
especially powerful one because it invokes well-known concepts. It 
simply applies the concepts of steady state, thermal equilibrium, and 
Fermi level to insulators and shows that the important features of 
space-charge layers in insulators can be described in terms of only one 
parameter, the work function or Fermi level position in the insulator 
at thermal equilibrium. (A more detailed discussion is available in 
the work of Simmons.*) With this as a starting point, it is possible to 
discuss a wide variety of charging phenomena in insulating thin films 
in an intuitively understandable way; and it should provide a basis 
for more quantitative analyses of a number of such effects. In the 
following section, the concepts discussed here are applied to the 
silicon /SiO2/Al,03/metal structure with emphasis on the space-charge 
region that builds up near the insulator/insulator interface during 
deposition of the Al:O3. 


IV. THE Si0./AI,0; SYSTEM MODEL 


The concepts of the preceding section can be applied to the double- 
layer structure shown in Fig. 2 in which the first layer is SiOz and the 
second is Al.O3. Previous experimental results have indicated that the 
density of trap levels in thermally grown SiO: is sufficiently low that the 
Debye length should be much larger than the typical SiO, film thickness 
of a few thousand Angstroms.’ This supports the assumption made 
previously that no finite charge density exists in the bulk of the SiO. 
film. On the other hand, the trap density in deposited Al,O; films has 
been found to be much larger so that it is reasonable to assume that 
the Debye length is small compared to the Al.O; film thicknesses that 
we will consider.§- The Al,O; films of interest are deposited at 900°C. 
The double-insulator system is assumed to be at an elevated tempera- 
ture during the deposition for a sufficient length of time that thermal 
equilibrium will be established, and under these conditions, the poten- 
tial diagram versus position will be like that shown in Fig. 5. Since the 
Fermi level in the silicon must coincide with that in the Al.O; at 
thermal equilibrium, a contact potential difference will exist that will 
result in (2) an electric field in the SiOe, (22) a space-charge region in 
the Si at the Si/SiO2 interface, (277) a space-charge region in the Al.O; 
at the SiO2/Al.O3 interface, and (iv) zero electric field in the bulk of 
the Al.O3; film. In addition, depending on the Al.O; surface boundary 
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Fig. 5—Band diagram of the Al.O;/SiO2/Si system in equilibrium at high tem- 
perature. 


conditions, a space-charge region may exist near the outer Al.O; 
surface. 

The important parameter in the SiO2/Al,03 system is the total 
potential difference V. that must build up to align the Fermi levels 
in the silicon and the Al.O;. This potential difference will be made up 
of a potential V.. across the SiO» and the drops in potential due to the 
band-bending regions in the silicon near the Si/Si0:2 interface and in the 
Al.O3 near the SiO02/Al.O; interface, AV, and AV 4, respectively (see 
Vig. 5). That is, 

Ve = Veo aa AV; + AV a, (8a) 


where the system equilibrium condition is given by 
(E¢/2) + oa + Ve — bi — ba = 0. (8b) 


If it is assumed that the work function argument applies to the 
8i02/Al.03 system, then 


Ve = Xa + ba — X; — (Ee /2), (9) 


where X4 is the AlzO3; electron affinity at the deposition temperature; 
ga iS the energy separation between the Al.O; conduction band and 
the Fermi level in the bulk; x; is the silicon electron affinity; and the 
Eg/2 term represents the assumption that the silicon is intrinsic at the 
elevated temperature of interest. 
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Since the electric field in the SiO. is determined directly by the 
potential drop V., and the oxide thickness 7',,, the net charge present 
in the narrow space-charge region at the Alz:03/SiO2 interface is 
given by 

Qii = (€02/T o2){Ve — AVs — AVa}. (10) 


Over a wide range of SiO» thickness, AV, and AV, will be negligible 
compared to V,. For example, even assuming that V,, corresponds to 
an electric field in the SiOz of 10° volts/em, AV, is less than 0.1 volt at 
900°C.* Thus, eq. (10) can be approximated by 


Qi = (€ce/ Toa) V. ~y en Pen, + PA = x = (H¢/2)}. (11) 


From eqs. (10) and (11), several interesting properties of Q;; are 
apparent. First, its magnitude is relatively independent of the quality 
and reproducibility of the Al,O3 provided only that ¢4 + X4 is repro- 
ducible and AV, is negligible. The band-bending AV4 may vary 
markedly from sample to sample depending on the trap density, but as 
long as the space-charge region is relatively narrow so that AV 4 is small 
compared to V., this variation will have no significant effect. Second, 
the magnitude of Q;; varies inversely with the SiOz thickness 7',. and 
is independent of Al,O; thickness T'4. This effect provides a straight- 
forward and unique prediction of the model that can easily be tested 
experimentally. 

If the system is now cooled to room temperature, the conductivity 
of the Al.O; will be reduced to the point where the space-charge 
regions will not move or change under application of an electric field 
for long periods of time, and these regions will be effectively frozen in. 
We must consider the two possible space-charge regions in the Al2.Os, 
one at the Al,03/SiO2 interface, and the other at the outer surface. 
The contributions V;; and V,,, respectively, of these charge layers to 
the flatband voltage is given by (see eq. 2) 


ToxtT —_ 
Vii = _ - | teers | pii(a)dx = Da it (12) 
Tor €A €A 
and 
ToztTA _— 
Wes | Fe Ea | oma)ae, (13) 
Tox €A 


where pii(%) and pn(x) are the net charge densities at the SiO2/A1.03 
and Al,0;/metal interfaces, respectively. In the limit where p;;(x) is 
located in a plane at the interface, an effective interface charge density 
Q:; is defined by eq. (12). The term V, is a constant independent 


*This can be shown from an integration of Poisson’s equation and using the fact 
that at 900°C the intrinsic charge density in the silicon is approximately 10 cm=°. 
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of T.. and 7'4 if pm(x) is a function only of the distance (7.2 + Ta — x) 
between the charge and the metal. Combining (2), (11), (12), and (18) 
gives for the flatband voltage of the double-insulator structure 


Vra(D.1.) = (¢m,a + bi: — bs) + Vin 
hs Oss Pow -+- we | = V. conta ) (14) 


€or €A €A T ox 








and, for completeness, the flatband voltage of the single insulator 
[eq. (3) ] is 
Qss 


€ox 


Vra(8.1.) _ (dm, oz = ds) _ 


where ¢, = $2 + 4p. 


fips (15) 


V. EXPERIMENTAL RESULTS 


5.1 Preliminary remarks 


Dual-dielectric mis capacitor structures with various insulator 
thicknesses were fabricated on n- and p-type silicon substrates with 
resistivities of approximately 10 ohm-cm. For the n-type substrates 
both (100) and (111) orientations were investigated. The SiO». was 
thermally grown at 1100°C using oxygen bubbled through 80°C water 
as the ambient. The Al.O3; was vapor deposited on the SiOz at 900°C 
from an AICl; source. The details of the AleO; deposition process are 
given in Ref. 1. The insulator thicknesses, 7',, and T'4, were varied by 
varying the growth and deposition times of the SiOz and the Al.Os, 
respectively. 

The deposition of the Al,O; duplicated exactly the procedure used 
for fabricating integrated circuits. As such, a layer of SiOe was also 
deposited on top of the Al,O3, which in the fabrication of integrated 
circuits is used as an etch mask for defining patterns in the Al.Q; 
film. For our samples, this layer of Si02 was chemically stripped prior 
to any measurements or any further processing. 

One feature that may be important in this study is the method of 
formation of the metal electrodes. Depending on the deposition 
technique used, the samples may be heated for a sufficient time during 
the metal deposition to form a space-charge region at the Al,O;/metal 
interface. However, if this induced space charge is reproducible and 
constant from sample to sample and is spatially constrained to a region 
very near the interface, it will only influence the flatband voltage via 
the constant voltage term V, in (14). Experimentally, we shall 
attempt to assure the reproducibility of this possible space-charge 
effect by measuring the mis structures at room temperature with a 
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mercury electrode." Several samples were also investigated with 
thermally evaporated titanium-aluminum electrodes. 

Initially, we attempted to vary the Al.0; thickness by etching in 
discrete steps rather than by varying the deposition time. This ap- 
proach was abandoned because of nonuniform etching of the Al.O3. 
In Fig. 6, scanning electron micrographs are given of the surface of 
the Al,O; as deposited and after etching a portion of the layer. Micro- 
scopic thickness variations (+500 A) are evident after etching. Since 
these variations lead to significant errors in the flatband voltage 
measurements, the Al,O3 thickness was varied only by varying the 
growth time. 

Experimental data for each sample investigated were obtained by 
means of high-frequency capacitance-voltage (C-—V) analysis.!2 Such 
measurements, obtained using either the mercury probe electrode or 
thermally evaporated titanium-aluminum thin-film electrodes, enable 
one to obtain accurate measurements of the insulator thickness and 










AS DEPOSITED 5 MIN ETCH-BACK 
Ta © 1750 A Ta © 1250 A 





12 MIN ETCH—BACK 17 MIN ETCH—BACK 
Ta = 500 A Ta OA 


Fig. 6—sem photographs depicting the increasing roughness of the Al.O; surface 
as it is etched back using phosphoric acid. 
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the associated flatband voltage. Each data point reported is the average 
flatband voltage calculated from at least three measurements per- 
formed on each sample. The typical spread in measurements is 0.10 volt. 
When measurements are performed on the double-insulator (Al,03/ 
SiO.) structure, a measurement of Vrg(D.I.) Leq. (14) ] is obtained. 
The normalized accumulation capacitance Cacc (farads/em?), which 
is a measure of the insulator thicknesses, is given by 

—l 

Cw(D.I) = | 2+ 24 |". 


€ox €A 


(16) 


When the Al,O3 is completely etched off and C-V analysis is con- 
ducted on the remaining single insulator (SiOz), then measurements 
of Vrp(S8.1.) [see eq. (15) ] are obtained. The normalized accumulation 
capacitance in this case yields a measurement of 7’. since 





€ox 
Cacc(S.1.) = 7 (17) 
By combining eqs. (16) and (17), accurate measurements of both T'4 
and 7'o2 are obtained. 

It is possible to obtain independent quantitative values for Q,, for 
each sample studied after the Al.O; is etched off if the constant term 
(bm,.or — $s) In eq. (14) is known. Measurement of (¢n,02 — $s) can be 
accomplished by SiOs etch-back experiments in which V,r,(8.I.) is 
measured as the SiO» layer is successively thinned by etching in a 
dilute hydrofluoric acid solution. Typical data obtained with a mercury 
probe on four different samples are shown in Fig. 7. As expected, there 
is a linear relationship between Vr,(S.I.) and 7'.z and an extrapolation 
of this relationship back to T., = 0 indicates that 0.67 volt is the 
appropriate value of (¢mn,oc — $s) for n-type, 10-ohm-cm Si and a 
mercury electrode. This value is in excellent agreement with previously 
determined values.1* The experiment also provides independent 
verification of the assumption that there is negligible space charge 
in the bulk of the SiOz. 

Based upon the value of Q,. for each sample, it is possible to charac- 
terize the flatband voltage shift due to the Al,O;. Correcting for the 
Q.; term and, additionally, subtracting the constant term (¢n,02 — $s) 
from eq. (14), a corrected differential flatband voltage AVrg can be 
defined as: 


AV rz 





+ 


Pgs F- 
Vra(D.1.) = Doves 5 os) + Oss é mal 


(dm,A + Pit =: dm, ox) + V. ( =e )( 7) + V s- (18) 


€A 
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Fig. 7—Plot of the Hg/Si02/Si flatband voltage for four wafers as a function of 
the Si02 thickness. Data obtained from SiO2 etch-back experiments on 10-ohm-cm, 
n-type, (100), Si substrates. 


5.2 Results pertaining to Qi: 


Experimental values of Vrs(D.I.) for the dual-insulator structure 
are plotted in Fig. 8 as a function of the Al.O3 thickness 74 for a SiO» 
thickness of 1200 A on n-type, (100) substrates. These results were 
obtained with a mercury electrode. Although there is considerable 
scatter in the data, it is clear that Vrg(D.I.) increases monotonically 
with increasing Al.O; thickness which is in agreement with eq. (14) if 
the sign of V. is such that a net negative charge exists at the Si02/Al.03 
interface. The experimental uncertainties in the Vrg measurements 
are estimated to be +0.05 volt. Correcting this data for Q,. and sub- 
tracting (¢@m,oc — $s), the results are replotted in Fig. 9. This refine- 
ment technique leads to a considerable reduction in the scatter in the 
data and demonstrates that AV rz is a linear function of the Al.O3 
thickness 7'4, as predicted by eq. (18). The linear relationship also 
provides striking evidence that Q::, the negative charge at the Al,O3/ 
SiO, interface, is constant from sample to sample for Al2O3 thicknesses 
in the range of 500 A to 2500 A if the SiO. thickness is held constant. 
This is in agreement with the postulated model and provides experi- 
mental verification of the assumption that the Debye length in Al,O3; 
is small compared to the Al.O; thickness. Given that the Debye length 
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Fig. 8—Plot of the Hg/Al,03/SiO2/Si flatband voltage as a function of the Al,O; 
thickness for a constant SiO thickness (7.2 = 1200 A) on n-type, (100), Si substrates. 


is small compared to 500 A at 900°C, estimates of AV, assuming 
charge densities in excess of 10'8 cm7? indicate that AV4 will be less 
than 0.1 volt and, thus, negligible as previously assumed. 

The data presented so far proves that Q;: is negative and a constant 
for fixed SiO. thickness independent of Al.O3 thickness. Another 
prediction of our model is that Q;; is inversely proportional to the SiO» 
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Fig. 9—Plot of the corrected differential flatband voltage as a function of the Al,O; 
thickness for a constant SiO thickness (7., = 1200 A) on n-type, (100), Sisubstrates. 
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Fig. 10—Plot of the corrected differential flatband voltage as a function of the 
Al.O; thickness for various values of SiO2 thickness on n-type, (100), Si substrates. 


thickness [eq. (11) _]. That this is indeed the case is shown by the data 


presented in Figs. 10 and 11, which were obtained with a mercury 
probe and are for (100), n-type, silicon substrates. In Fig. 10, AV rz is 
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Fig. 11—Plot of the corrected differential flatband voltage as a function of the ratio 
of Al.O; thickness to SiO2 thickness on n-type, (100), Si substrates. 
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Fig. 12—Plot of the corrected differential flatband voltage as a function of the 
ratio of AleO3 thickness to SiOz thickness on n-type, (111), Si substrates. 


plotted versus Al,O; thickness for three different SiO. thicknesses. In 
each case, a linear relationship between AV rz and Al,O; thickness is 
obtained, and the increased slope obtained with smaller SiO, thick- 
nesses indicates that Q,; does increase as the SiO» thickness is decreased. 
The data of Fig. 10 are replotted in Fig. 11 as a function of T4/T oz, 
the ratio of AlsO; thickness to SiO. thickness. As expected from eq. 
(18), the AVrg versus T74/T 2 relationship is accurately represented 
by a straight line over the 74/T.2: range investigated (0.5 to 8), 
indicating that Q;; is inversely proportional to the SiO» thickness. The 
slope of the straight line in Fig. 11 corresponds to a V, value of 0.88 
volt.” This value for V, was obtained in all the measurements on (100) 
substrates within +0.1 volt. 

Similar measurements were also performed with (111) oriented, 
n-type substrates. The larger values of Q,, inherent in the (111) 
orientation meant that the Q,, correction factor was much larger and, 
hence, the accuracy of the results was somewhat poorer. Results for 
(111) samples are given in Fig. 12, where AV rz is plotted as a function 
of the T4/T 2 ratio. Again these data were obtained with a mercury 
probe. The straight line shown in Fig. 12 is a best fit to the data if 
the slope of the line is restricted to correspond to a V, value of 0.88 volt. 
Considering the possible errors due to the Q,, correction, the straight- 
line fit of the data in Fig. 12 is good enough to conclude that the 
value of Q:; is independent of the substrate orientation for the two 
orientations investigated, (100) and (111), and in complete agreement 
with the predictions of our model. 


* For 1000 A of SiOz, a V, value of 0.88 volt corresponds to a Q;; value of 1.9 & 10" 
charges/cm?. 
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5.3 Results pertaining to Q.. 


All of the experimental results presented so far have focused on 
Q::, the charge at the Al,O3/SiO2 interface. Some interesting facets 
of Q;s, the charge at the SiO;/Si interface, were also discovered during 
our study and are discussed in the following. As mentioned previously, 
a value of Q;,; was determined for all samples by measuring the flatband 
voltage Vrz(S.1.) of the single-insulator structure after removing the 
Al,O3 and then calculating Q,, using the (@m,oz — $s) value for the Hg/ 
SiO2/Si system as determined in Fig. 7. A plot of Vrg(8.I.) versus 
T or, the SiOz thickness, for n-type, (100) substrates is given in Fig. 13. 
Previous results published in the literature have shown that for single- 
insulator (Si02/8i) structures, Q,; 1s independent of the SiOz thickness."4 
If this were the case for our structures, we would expect to find a 
linear relationship between Vrz(S.I.) and 7... with an intercept on the 
Vra(8.1.) axis equal to (¢m.or — ¢s) = 0.67 volt. The results given 
in Fig. 13 indicate that this is not the case. Although the data could 
be interpreted as being consistent with a linear relationship, they are 
definitely not consistent with an intercept equal to 0.67 volt. The results 
are more consistent with the supposition that, to first order, Vrz(8.1.) 
is independent of T oz. 

A more detailed study of Q,; was pursued by preparing samples of 
various SiOz thicknesses (n-type, Si, (100)) and measuring Vr,(8.1.) 
for each sample. Approximately 500 A of Al;O3; was then deposited 
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Fig. 13—Plot of the Hg probe single-insulator flatband voltage after AlxO; deposi- 
tion and etch-off as a function of SiO2 thickness for n-type, (100), Si substrates. 


DOUBLE-INSULATOR MOS 705 


on all samples, the AlzO; was removed by etching, and Vr ,(S8.I.) was 
remeasured. Finally, for each of the samples, the SiO. was etched 
back in steps, and Vrz(8.I.) was determined as a function of SiO: 
thickness. The results are given in Fig. 14. 

Prior to Al,O; deposition, the results are consistent with the samples 
having a constant value of Q,, independent of 7.2, that is, a linear 
relationship exists between Vr,(S8.1.) and T.. with an intercept equal 
to 0.67 volt. However, after Al,O; deposition, Vrz(8.I1.) is seen to be 
essentially independent of 7’... Furthermore, if the SiO2 is now etched 
back, a linear relationship between Vyg(8.I.) and SiO, thickness is 
obtained with an intercept equal to 0.67 volt. In Fig. 14, results of the 
etch-back study are given for only one representative sample, since the 
results obtained on the other samples were similar. 

The conclusion which follows from the results given in Fig. 14 
is that before the deposition of the Al.03, Q,, is independent of T'o:, 
whereas after deposition, the value of Q,, is changed, the amount of 
change depending upon 7’,,, the SiO» thickness. This effect is further 
illustrated by the results given in Fig. 15, where Q,, after AleO3 deposi- 
tion and etch-off is plotted versus 1/7’... for both (100) and (111), 
n-type substrates. For both orientations, the data are seen to fall 
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Fig. 14—Plot of the single-insulator flatband voltage before and after Al.O; 
deposition and after final SiOz etch-back, indicating the change in Q,; induced by the 
Al.O; deposition on n-type, (100), Si substrates. 
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Fig. 15—Plot of Q;s after Al2O; deposition as a function of the SiO» thickness for 
n-type, (100) and (111), Si substrates. 


along a straight line consistent with the equation 
Oss — (Vo€ox/ T oz) + Giro: (19) 


The slopes of the two lines in Fig. 15 were taken to be equal to each 
other (v. = 0.38 volt), and it is observed that an excellent fit to the 
two sets of data is obtained with the one v, value. The background 
charge densities, Q,;., or equivalently the values of Q,, for large values 
of T.2 are approximately 1.0 X 10" and 1.0 X 10 charges/cm? for 
(111) and (100) orientations, peboenvely: For the (100) orientation, 
Qsso 18 negligible. 

The results presented so far have established that for the double- 
insulator structure, both Q;:; and Q,, depend on 7',, and that these 
dependencies can be written as: 


Qi: = oe Ve/ Taz 
Oss 7 Qsso = Cardo) Tony 


where V, ~ 0.88 volt and v, ~ 0.88 volt. Thus, the single-insulator 
flatband voltage after Al,O; deposition can be written as 


Vra(8.1.) zs Ca Ss ds) — Vo — CT cat on) Devos (21) 


(20) 
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which is consistent with the after-deposition results presented in 
Fig. 14. Considering Q,,. as a fundamental property of the SiO./Si 
interface, which is not affected by the Al,O; deposition, it follows that 
the contribution to the single-insulator flatband voltage, 6V rz, induced 
by the portion of Q,, that is influenced by the Al,0; deposition is 


OVrg = UV. (22) 


The voltage drop across the SiO» during the deposition of the Al,O3 
is V., and if this is considered as a stress voltage applied to the Si02/Si 
interface, then 

V.(stress)/6Vrp =a = 2.3. (23) 


It is interesting to note that this result is in good agreement with 
previously published results relating stress voltage to saturated 
flatband voltage shift for SiO2/Si structures. Specifically, in Ref. 15 
it was observed that the ratio of stress voltage to saturated flatband 
voltage shift was given by: 


3.33 at 350°C 


ag bee at 450°C. 24) 


One implication of the relationship given in eq. 20 for Q,; is that for 
(100) substrates, where Q,;. is negligible, the flatband voltage of a 
double-insulator structure will be insensitive to a SiO2 thickness 
variation. Thus, in MOSFET integrated circuits with a (100) substrate, 
where a thick SiO, layer is used to inhibit parasitic inversion, the 
contribution of the Q,, term to the parasitic inversion voltage will be 
independent of SiOz thickness. 


5.4 Results pertaining to p-type substrates 


The substrate conductivity type does not appear directly in the 
model that has been proposed for the magnitude and origin of Q;:, and 
all of the results presented so far have been for n-type substrates. Since 
the silicon substrate will be intrinsic at 900°C, the deposition tempera- 
ture of the Al.O;3, the voltage drop across the oxide V, will be the same 
for both n-type and p-type substrates and, thus, Q;; at 900°C should 
also be independent of the conductivity type of the substrate. If, 
during the cool down after Al.O; deposition, Q;; is frozen in at a tem- 
perature at which the silicon is still intrinsic, then the value of Qi: 
measured at room temperature should not depend on whether the 
substrate is n-type or p-type. Results obtained with (100), p-type 
substrates are presented below. 

A plot for p-type substrates similar to that of Fig. 11 (for n-type 
substrates) is given in Fig. 16, where AV rz is plotted as a function of 
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uw 


AVrg IN VOLTS 


Hg PROBE 
p—TYPE Si<100> 





RATIO (Ta/Tox) 


Fig. 16—Plot of the corrected differentia] flatband voltage as a function of the ratio 
of Al,O; thickness to SiOz thickness on p-type, (100), Si substrates. The solid line is 
taken from Fig. 11. 


T/T. The solid line corresponds to the linear fit of the data of 
Fig. 11. Although there is general agreement between the results for 
the n-type substrate (solid line) and this experimental data, the large 
amount of scatter in the data must be recognized. Figure 17 is a plot 
of the single-insulator flatband voltage for the p-type substrates after 
the AlsO; etch-off. The (¢m,o2 — ¢;) value was obtained by etch-back 
experiments, as outlined previously. Here again, a large amount of 
scatter in the data is evident. 

The variations in the above sets of data are not random scatter, 
but are due to some mechanism unique to the p-type substrates. In 
Fig. 18, the single-insulator flatband voltage Vrx(8.I.) is plotted as a 
function of 74 for n-type samples after Al2O; etch-off, and it is clear 
that no dependence on 74 or Tz is evident. In Fig. 19, similar data 
are plotted for the p-type samples and it is evident that in this case 
there is a dependence on 74, but again, no dependence on 7',;. For 
both cases, the conclusions regarding 7',, are obtained from the points 
in Figs. 18 and 19, which explicitly denote points of constant T'..2. 
Although a detailed explanation of this effect cannot be given, it 
is felt that this effect is due to the fact that boron-doped p-type 
wafers were used in the experiment. It is known that boron will greatly 
out-diffuse from a silicon substrate into an $102 layer in the presence 
of a high-temperature, hydrogen-containing ambient.!*17 Additionally, 
the introduction of this impurity into the SiOz may enhance its conduc- 
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Fig. 17—Plot of the single-insulator flatband voltage after AlO; deposition and 


etch-off as a function of the SiO: thickness for p-type, (100), Si substrates. 
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Fig. 18—Plot of the single-insulator flatband voltage after AlxO; deposition and 


etch-off as a function of the Al,0; thickness for n-type, (100), Si substrates. 
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Fig. 19—Plot of the single-insulator flatband voltage after Al,O3; deposition and 
etch-off as a function of the Al,O; thickness for p-type, (100), Si substrates. 


tivity at high temperature. The net result will be a drop in the effective 
stress voltage V, across the SiO, layer as a function of Al,O; deposition 
time* and a lowering of the Q;; and Q,s term. 

This hypothesis is consistent with the following data. In Fig. 20, 
distributions of the potential drop v, due to the induced Q,, term (see 
eqs. 19 and 20) are plotted for both the n-type and p-type (100) 
samples. It is observed that 


(t) Lower voltage drops for p-type samples occur than observed 
for the n-type samples (the lower values are correlated to the 
thicker Al,O; deposition). 

(zi) No zero (or negative) voltage drops occur. 

(227) The upper range of v, for the p-type samples (the lowest SiO. 
conductivity region), are bounded by the v, values observed 
for the n-type samples. 


One final point can also be made to support the hypothesis. It is 
possible to calculate Q;; values from the experimental data (via eqs. 
10 and 18) if an intercept voltage value (i.e., 774 = 0) is assumed; and 
from the experimental data for n-type samples, an intercept value of 


* In all experiments, the Al,O; deposition rate was a constant (75 A/min.). 
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Fig. 20—Distribution plot of the single-insulator flatband voltage observed on 
n-type and p-type, (100), Si substrates after Al2O; deposition and etch-off. 


1.25 volts was obtained. Additionally, Q,; values for each sample may be 
calculated. If the above hypothesis concerning a lowering of V, is 
correct, then the ratio of Q:; and Q,, (given in eqs. 20 and 23) should 
be independent of the absolute magnitude of V, for (100) substrates 
in which Q,,. is negligible. That is, 


Qe — Coit Toe = €or Ve /T oz 
and (25) 
Qii/ Qss = a. 

Figure 21 is a plot of Q:: vs Q.. for both the n-type and p-type samples. 
It is noted that a linear relationship exists for both sets of data with a 
the same for both (a 2.5). The several data points that deviate 
from the linear relation are associated with very small Q,. and Qi: 
values, and the deviation is most likely due to small errors in flatband 
voltage measurements. 


5.5 Results pertaining to observed potential jumps 


Extrapolating the linear relationships in Figs. 9, 11, and 12 to 
Ts = 0 gives a value: 


(dma + Git — dmoz) + Vn & AG + Vn 1.25 volts. (26) 
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While the experiments cannot determine the relative contributions of 
V., and A@ to the intercept value, it is worth pointing out the conse- 
quences of the two limiting possibilities. First, if V, is zero, then Ad 
is non-zero and equal to 1.25 volts. This implies that the work function 
model for barrier heights must be incorrect (otherwise Ad = 0). The 
second limiting possibility is that A¢é = 0 and V,, is non-zero. In this 
case, there must be a 1.25-volt band-bending effect at the outer Al,O3 
interface. Although we have not been able to perform an experiment 
that unequivocally separates the contributions of A@ and V,, to the 
intercept value, it is worthwhile to consider some additional items 
of relevant experimental information. 

First, it was stated previously that all samples are fabricated with a 
deposited SiO» layer on top of the Al.O3 layer. An etch-back experi- 
ment was conducted on the deposited SiO» (using n-type Si, (100), 
To: = 600 A, 74 = 500 A) and the flatband voltage was measured as a 
function of the equivalent SiO» thickness 7'.,: 


Ts => Legs + (€o2/€4)T a + T sio., (27) 
where 7'sio, equals the deposited SiO» thickness. The results of this 


1012 
8 


oO 


eS 






nN 


= 
(o) 
a 


@= n—TYPE <100> 
A= p— TYPE <100> 


CHARGE AT Al,03—SiO. INTERFACE, Q;, (CHARGES /cm2) 


2 
1910 2 4 6 8 qQt 2 4 6 8 4912 
CHARGE AT Si0—Si INTERFACE, Ogg(QUANTITY/cm2?)} 


Fig. 21—Plot of Q,:, the charge at the Al.0;/SiOe interface, versus Qs, the charge 
at the SiO2/Si interface, for both n-type and p-type, (100), Si substrates. 
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experiment are plotted in Fig. 22. It is evident from the linearity of the 
flatband voltage that no significant charge density is present in the 
bulk of this deposited SiOz. It is interesting to note that a positive 
potential jump of 1.22 volts in the flatband voltage is associated with 
the outer Al.O3/Si02 deposited interface (using the Hg metal electrode). 
This value is close to the 1.25-volt potential jump associated with the 
inner Al.O3/SiO.2 (thermal) interface. 
It is straightforward to show that the potential jump observed in 
Fig. 22 is given by 
PJ = (6m,4 + biti — bm,oz) (28) 


if it is assumed that there is no change in the charge distribution in the 
Al.O3 film when the deposited SiOz is completely removed, and that 
the barrier height of metal-to-deposited-SiO, is the same as the barrier 
height of metal-to-thermal-SiO2. With these assumptions, the con- 
clusion follows that V,, 0 and A¢ ~ 1.25 volts. 

Second, measurements were also made with titanium-aluminum 
evaporated electrodes. The double-insulator flatband voltage for this 
- metallization system is plotted in Fig. 23 as a function of the Al.O3 
thickness 74 for a constant SiO. thickness 7.2 1200 A (n-type Si, 
(100)). It is evident that the scatter in the data is much greater than 
that found for the Hg metallization. Attempts to refine the data proved 


n-TYPE <100> Si 
102 —cm 


Hg PROBE 
¢ PJ = 1.22+0.07 V 
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0 1 2 3 4 5 6, 7 8 
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Fig. 22—Plot of the triple-insulator flatband voltage, obtained by means of an 
etch-back experiment, as a function of the equivalent insulator thickness. 
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Fig. 23—Plot of the titanium-gate double-insulator flatband voltage as a function 
of the Al.O; thickness for a constant SiOz thickness (7.2 = 1200 A) on n-type, 
(100), Si substrates. 


fruitless due to the additional scatter observed in the single-insulator 
(SiO.2) flatband voltages (see Fig. 24). 

Some general comments can be made, however, concerning these 
data. The scatter in the double-insulator flatband voltage (Fig. 23) 
decreases with descreasing Al,O0; thickness, indicating that an un- 
controlled charging effect must take place in the Al.O; during the 
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Fig. 24—Distribution plot of the titanium-gate single-insulator flatband voltage 
after Al,O3 etch-off as observed on the samples reported in Fig. 23. 
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metal deposition. Extrapolating to T4 = 0 implies 
0.40 V < [Ves(D.1.)|rs-0] < 0.60 volt. (29) 


The distribution in the single-insulator flatband voltage Vrg(8.I.) can 
be characterized by (see Fig. 24) 


Average Vr,(8.1.) = — 0.48 volt 
Standard Deviation Vrz(8.I.) = 0.10 volt. 


Thus, a positive potential jump of 0.98 + 0.14 volt can be associated 
with the addition of the Al,O; layer. This value is in reasonable agree- 
ment with the previous potential shift results found for the Hg metalli- 
zation system. 

Third, in Ref. 4, the authors find on a Vrg(D.I.) versus 7'4 plot for 
various JT, values an intercept value of —0.80 volt at 74 = 0. This 
result, obtained on 2-ohm-cm, p-type, (100), Si substrates, is inter- 
preted by the authors to be the expected metal-to-silicon work function 
difference when using aluminum electrodes, a conclusion that may not 
be valid since Q,; measurements after Al,.O3 deposition were not 
reported. We shall now show that the results in Ref. 4 are in excellent 
agreement with our results by taking the results we have obtained 
with a Hg electrode and converting them to the results we would have 
obtained if we had used an Al electrode. 

To explicitly denote the use of a Hg electrode, eq. (26) is rewritten as 


(ue, A + ii => ug, oz) + Vix = 1.25 V (30) 
and from our measurements on p-type material 
ug, oz — %s = 0.08 V. (31) 


The difference in barrier heights between Hg and Al on the type of 
A],03 studied in this paper has been reported by Nigh!® and is given by 


ouz,aA — Gaia = 17 V. (32) 
Combining eqs. (380), (81), and (82) yields 
(pai,a + bis — be) + Vm = — 0.87 V. (33) 
The intercept value for T'4 = 0 predicted by (14) is given by 
Intercept (T4 = 0) = (¢a1,a + $i — bs) + Vn — Gale (34) 


Using the values », = 0.38 volt and Q,;. = 0 for (100) material in eq. 
(20), and combining eqs. (20), (83), and (34) yields 


Intercept (74 = 0) = — 0.75 V. 
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Thus, our corresponding intercept value for Al electrodes is almost 
identical to that reported in Ref. 4, and it strongly implies that the 
electrical properties of the different Al,O3/S10. films which determine 
the flatband voltage in mos structures are identical. More specifically, 
it indicates that the Q,, dependence on 7',, reported in this paper is 
also true for the structures studied in Ref. 4, and that the value of 
A@d + Vin both cases is the same. 


VI. SUMMARY AND DISCUSSION 


By varying the thicknesses of both insulators in a silicon/SiO2/ 
Al,O3/mercury Mos structure and accurately measuring changes in 
flatband voltage, we have established that a net negative space charge 
exists near the Si0./A1.03 interface, which is spatially constrained to a 
region much less than 500 A thick. The magnitude of this negative 
charge varies inversely with SiO. thickness and is the same for both 
(100) and (111) oriented n-type and p-type silicon substrates. These 
results are consistent with a model for space-charge formation based 
on work by Simmons on metal-insulator-metal structures. At the 
elevated temperature (900°C) of Al.O; deposition, the Al,O; is a 
good enough conductor that thermal equilibrium is established. 
Since the electrostatic screening or Debye length in Al,O;3 at this 
temperature is small compared to the Al,O; thickness of interest, the 
bulk of the Al.O; is at zero electric field, and a Fermi level can be 
defined that must align with the Fermi level in the silicon substrate. 
This requires that a fixed “contact potential”, experimentally found 
to be 0.88 volt, must exist across the SiOz at 900°C. The electric field 
associated with this potential generates a net negative space-charge 
layer near the SiO2/Al.03 interface. When the structure is cooled to 
room temperature, the conductivity of the Al.O3 reduces to a negligible 
value and the space charge is frozen in. The net negative charge can 
thus be considered to be the charge on the SiO: capacitance associated 
with the constant 0.88 volt contact potential. 

When the double-insulator flatband voltage is corrected for the 
independently measured Si/SiO» interface charge Q;, and the barrier 
heights of the single-insulator system, a corrected differential flatband 
voltage is generated. Extrapolation of this function to zero Al,O3. 
thickness reveals a potential jump of approximately 1.25 volts when 
using a mercury electrode. Similarly, a potential jump of approxi- 
mately 1.0 volt is found with titanium-aluminum electrodes. The inter- 
facial barrier energies that contribute to these jumps are shown not 
to be derivable from a simple work function argument. 

The large amount of scatter observed in the data where titanium- 
aluminum electrodes are used, compared to the very consistent data 
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obtained with the mercury electrode, implies that thermal evaporation 
of a metal onto an Al,O; film introduces significant variation in the 
flatband voltage. This effect is most probably due to a charging 
phenomena that occurs during the transient heating of the sample 
during evaporation. 

Measurements of Q,; made on samples before and after Al.O; 
deposition revealed that during this deposition, the value of Q,. was 
changed. After Al,O; deposition, it was found that Q,, could be written 
as the sum of two terms. One term was a constant background charge 
density that was independent of SiO. thickness and that had the 
values of 1.0 X 10" and 1.0 X 10" charges/cm? for (111) and (100) 
oriented substrates, respectively. The other term was orientation 
independent and inversely proportional to the SiO. thickness, indi- 
cating that it is derivable from a constant contact potential that was 
experimentally determined to be 0.38 volt. Thus, the electric field that 
exists in the SiO» during the deposition of the Al,O; determines not 
only Q;: but also a portion of Qes. 

The proposed charging model was also found to be correct for p-type 
substrates except that an additional effect was uncovered in that the 
effective contact potential decreased with increasing Al,O; thickness. 
This effect may be due to boron penetration of the thermal SiOz layer 
and an associated increased electrical conductivity at high temperature. 

According to the model presented, a contact potential exists across 
the SiO» at the Al,O; deposition temperature, which results in an 
electric field in a direction to drive mobile positive ions away from the 
$i/SiO, interface. This means that if some mechanism exists for either 
removing or immobilizing these positive ions when they reach the 
$i02/Al,03 interface, the AleO3; deposition is expected to stabilize 
the MIs system against ionic drifts. Such a mechanism may indeed 
be present since HCl, a by-product of the Al,O; formation reaction, is 
known to be an excellent sodium getter. While the importance of this 
electric field in accounting for the stability of the 8i02/Al.03 system 
is not presently clear, it seems reasonable to assume that net positive 
charge at the insulator-insulator interface would make it much more 
difficult to remove or immobilize positive ions in the SiOe during 
second insulator deposition, since the positive ions would then tend 
to drift to the Si/SiOz2 interface. 

Since the presented model for space-charge layer formation at 
insulator-insulator interfaces is relatively insensitive to the nature of 
the deposited insulator, provided the assumption of thermal equi- 
librium at the deposition temperature is correct, considerations 
similar to those given in this paper for the 8i02/Al,0; system should 
apply to other SiO2/deposited insulator systems. 
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A 1-Watt, 6-Gigahertz IMPATT Amplifier for 
Short-Haul Radio Applications 


By J. E. MORRIS and J. W. GEWARTOWSKI 
(Manuscript received September 14, 1973) 


A 1-watt IMPATT diode amplifier has been developed for short-haul 
FM radio relay applications in the 6-GHz common-carrier band. The 
amplifier is used in the new TM-2 system and as part of a retrofit package 
to upgrade the performance of the existing TM-1 system. Amplification is 
provided by a single silicon IMPATT diode which is used in an injection- 
locked mode. A finned heat sink provides IMPATT diode cooling by 
natural air convection within the radio bay. The diode is expected to have 
a mean life greater than 10 years, and wt can be replaced in the field without 
the use of special tools or equipment. This microwave-integrated amplifier 
contains the rf samplers and detectors necessary to monitor both input 
and output rf power levels. The input power monitor also provides an 
input to a power-supply squelch circuit that removes dc power from the 
IMPATT diode tf the rf input signal level becomes too low for adequate 
performance. The influence of the system requirements upon the amplifier 
design 1s described, and data on system performance are presented. 


|. INTRODUCTION 


The mmpatr diode has been developed to the point where several 
watts of cw power can be generated reliably in the microwave fre- 
quency range. This negative-resistance device used in conjunction with 
a circulator comprises a reflection amplifier suitable as the power 
amplifier in a microwave communications transmitter. In the present 
application, the diode operates in the injection-locked oscillator mode. 
It was demonstrated by Tatsuguchi, Dietrich, and Swan that such an 
amplifier using a single silicon mmpatrT diode could meet the basic 
performance objectives of a typical short-haul radio-relay system.! 
The amplifier operates with a nominal gain of 20 dB and a noise figure 
of less than 52 dB. The corresponding system performance is better 
than 22 dBrncO per hop for a 1200-circuit message load. The amplifier’s 
system performance is found to be dominated by thermal noise, with 
intermodulation distortion negligible. The dc-to-rf efficiency is 4 
percent. 
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To be useful to the system, the amplifier package also contains rf 
samplers and detectors necessary to monitor the rf input and output 
power levels. The input power-monitor circuit furnishes the input 
information for a power-supply squelch circuit. If the input rf level 
drops low enough so that the locking bandwidth of the amplifier 
becomes small, the power supply is turned off, preventing the mpaTr 
oscillator from free-running out of the assigned frequency range. The 
de power is automatically restored when the input rf level returns to 
normal. The amplifier also contains harmonic suppression filters to 
prevent radiation of spurious tones. The amplifier has standard 
WR-159 waveguide input and output ports with vswrs of less than 
1.07 across the band. 

To be suitable for manufacture, an economical design was evolved 
based on the microwave integrated-circuit techniques successfully 
employed in the TH-3 system by Dietrich.? The construction consists of 
a thin-film strip-line pattern on a suspended alumina substrate, which 
is mounted in a die-cast aluminum housing, connected to a coaxial 
section containing the impatTr diode, the tuning mechanism, and a 
second harmonic filter. In addition, a wide range of tunability had to be 
incorporated to accommodate a wide range of diode parameters, both 
for initial manufacture and field replacement of the diode. Both 
frequency and output power adjustments are provided. All these 
features have been successfully accomplished in the amplifier designed 
for manufacture. 


li. AMPLIFIER DESIGN 


The amplifier design is based upon the use of a single silicon IMPATT 
diode used in a phase-locked oscillator mode. This mode of operation, 
described below, is chosen since it permits the relatively high gain of 
approximately 20 dB to be obtained stably in a single stage. 


2.1 Operating point selection 


The choice of an operating point for the amplifier follows the method 
described by Tatsuguchi et al.! Figure 1 illustrates typical contours 
of constant system thermal noise performance, in dBrncO per hop, 
plotted on coordinates of amplifier output power versus amplifier 
noise figure. The system performance contours shown apply to the 
highest frequency message slot of one particular short-haul rm system 
configuration that is operated at a 1200-message circuit loading. The 
contours assume that a +10-dBm level signal is available to drive 
the amplifier. This input power level is the minimum value anticipated 
in one of the systems in which this amplifier will be used. The options 
open to the amplifier circuit designer are illustrated on the same figure 
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Fig. 1—Contours of constant system thermal-noise performance, in dBrncO per 
hop, plotted on coordinates of amplifier output power versus amplifier noise figure. A 
particular amplifier’s performance is indicated by the solid lines at several imparT 
diode de currents. Typical performance obtained on a large sample of amplifiers when 
adjusted for 1-watt output power is shown. 


by the superimposed contours of mmpatr amplifier performance at 
various de power levels. For a given de power level, the operating 
point is a function of the microwave circuit impedance seen by the 
IMPATT device. The shape of these curves is due to the fact that an 
IMPATT device becomes noisier as the rf level is increased. From such 
curves, it becomes apparent that operation at the maximum possible 
rf power will result in poor system performance. Optimum performance 
occurs at neither maximum rf power nor minimum noise. It is instruc- 
tive to note that the optimum performance, i.e., lowest dBrncO 
number, occurs with the largest dc power. The use of high de powers 
must be tempered by reliability considerations, which generally dictate 
the use of lower powers. 

For this amplifier application, the trade-off between rf output power, 
FM noise, and diode reliability formed the basis of the decision to 
operate at 1-watt rf output with 24 watts of de supplied to the ImpaTT 
diode. At this operating point, the diode junction temperature is 
expected to be approximately 200°C in convection-cooled radio bays 
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operating in room ambient temperature up to 50°C. This operating 
point is expected to provide a mean diode life greater than 10 years. 
This reliability is the result of careful device processing combined with 
low thermal impedances both within the diode package and between 
the diode case and ambient air. 


2.2 Oscillator mode 


The impatr diode is operated in an injection-locked (phase-locked) 
oscillator mode, shown schematically in Fig. 2. The mpatt device and 
its associated resonating circuitry terminate one port of a circulator 
in a negative impedance. In the absence of an input signal, a free- 
running oscillation at frequency fy) occurs, which is coupled to the 
output through the circulator. When an appropriate input signal is 
added, the oscillation frequency locks to the input over a band of 
frequencies 2Af, approximately symmetrical about f». The free- 
running frequency is adjusted to the desired operating channel. Figure 
2 illustrates the power and phase variations that occur across the 
locking frequency band. The power levels and oscillator external 
Q (Q-x) are chosen such that Af is at least 10 times the highest modu- 
lating frequency. In this way, only the center linear portion of the 
phase variation curve is used, and phase distortion is minimized. 

Since the rf output power is fixed from other considerations (system 
performance and fading margin) and the rf input power available in 
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Fig. 2—Simplified representation of an injection-locked oscillator. 
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the systems in which the amplifier is to be used is limited, the designer 
is required to provide a circuit of the lowest possible Q. For this 
amplifier, a low Q circuit is provided by the use of a diode circuit 
consisting of a coaxial one-quarter-wavelength transformer plus a 
short section of coaxial line between the transformer and the diode 
that series-tunes the impaTT diode’s capacitance. With this resonator, 
the circuit Q is sufficiently low that Q.. is largely determined by the 
IMPATT device itself. 


2.3 IMPATT characteristics 


The rmpart diode used for this amplifier is an n-type silicon diode 
whose junction side is bonded to a metallized diamond within a copper 
and ceramic microwave pill package.* The large-signal rf charac- 
teristics of the diodes are measured near 6 GHz using the method 
described by Decker et al.t The diode wafer admittance is measured on 
all devices at 24-watts de, with an rf voltage corresponding to the 
diode’s operating point in the amplifier. Wafer susceptances are speci- 
fied at 19.0 millimhos. Tuning is provided to accommodate a range of 
diode susceptances. The wafer Q, defined as the magnitude of the ratio 
of wafer susceptance to wafer conductance, has values that vary by a 
factor of 2.5 to 1. 


2.4 Circuit description 


The requirements of practical radio-relay equipment dictate an 
amplifier circuit somewhat more complex than the simple circulator, 
diode, and resonator shown in Fig. 2. A more complete schematic of the 
amplifier is shown in Fig. 3. The circuit contains three circulators, 
of which the center circulator corresponds to the one shown in Fig. 2. 
Additional circulators with one port resistively terminated are used 
at both the amplifier’s input and output to provide isolation from the 
effects of external reflections and to provide input and output return 
losses better than 30 dB. 

The de power for the rmmpatTr diode from the current-regulated 
power supply is coupled to the oscillator port of the center circulator 
through a resistor and a band-stop filter tuned to 6 GHz. The resistor 
is used here to provide the high resistive impedance at low frequencies 
that has been shown by Brackett to prevent spurious oscillations.® 
The de power is isolated from the remaining rf circuit by a series 
capacitor in the main rf circuit adjacent to the band-stop bias filter. 

On the output side of the center circulator, a small sample of the 
amplified output is picked off by a nondirectional coupling probe. 
This sample of the rf output is detected using a point contact diode 
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Fig. 3—Schematic diagram of complete impaTrT diode amplifier. 


to provide a de current for the radio bay meter panel and alarm 
functions. 

On the input side of the center circulator, a sample of the input 
power is taken and detected for bay panel metering and to operate a 
power-supply squelch circuit. The power-supply squelch operates to 
remove dc power from the impart diode if the rf input to the amplifier 
falls below a prescribed level. This effectively prevents free-running 
oscillations by the mmpart oscillator, whose free-running frequency is 
not sufficiently stabilized to prevent interchannel interference. A 
directional coupler is used for the input coupler to provide a good 
match on the circulator common-arm and to provide, via its 20-dB 
directivity, discrimination against leakage of power generated by the 
IMPATT diode. 

A band-stop filter is used on the oscillator port of the center circu- 
lator to prevent second-harmonic energy generated by the impaTr 
from interfering with the operation of the monitor circuits. A low-pass 
filter is located at the amplifier’s output to ensure that all harmonics 
are suppressed. 


Ill. CIRCUIT FABRICATION AND TUNING 


Most of the circuit is fabricated using the microwave integrated- 
circuit techniques developed for use in the TH-3 system.? Film inte- 
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grated circuits (Fics) consisting of patterns defined photolithographi- 
cally on 0.024-inch (0.61-mm) thick unglazed alumina are used as a 
suspended-substrate strip-line transmission-line medium. The strip- 
line circuitry, as well as the amplifier’s waveguide input and output, 
are contained in a die-cast aluminum housing, shown in the amplifier 
photograph, Fig. 4. The mrparr diode and its resonator are contained 
in a short section of coaxial line that projects perpendicularly from the 
housing and is topped by the large, finned heat sink used for mpaTT 
diode cooling. Adjustments are provided on the coaxial section for 
field tuning of frequency and power output. 


3.1 Strip-line circuits 


The layout of the circuitry within the die-cast housing is illustrated 
in Fig. 5 and shown pictorially in Fig. 6. The ceramic substrates are 
located within a narrow channel to avoid multimoding problems. The 
complex substrate shape is fabricated by an automated laser-cutting 





Fig. 4—Complete amplifier. 
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Fig. 5—Arrangement of waveguide and strip-line circuit within die-cast housing. 


technique. The amplifier’s full-height waveguide input is shown on 
the right. A thin-film probe transition from the input waveguide to 
the suspended-substrate strip-line couples the input signal to the first 
of the three circulators. This circulator, with one port terminated in a 
thin-film resistor, provides the necessary input isolation. The circulator 
and termination designs follow closely those described by Dietrich,? 
modified to improve the temperature stability. A directional coupler, 
located between the input and center circulators, diverts approximately 
10 percent of the input rf signal to the input detector diode to generate 
the de needed for bay panel metering and power-supply squelch 
functions. 

The remaining input signal is coupled to the oscillator port of the 
center circulator. The series capacitor, which dc-isolates this port, is 
realized by a narrow, meandcring, interdigital gap in the thin-film 
conductor. The 6-GHz band-stop bias filter is realized by a high- 
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characteristic-impedance line along which, at quarter-wavelength 
intervals, are placed two quarter-wavelength-long open-circuited stubs 
of lower characteristic impedance. The bias circuit is completed by 
a 47-ohm power resistor that is clamped to the aluminum housing to 
maximize heat transfer. Several ferrite beads are placed on the leads 
of this resistor to provide additional stability against bias circuit 
oscillations. 

An open-circuited stub, one-half-wavelength long at 6 GHz, which 
connects to the oscillator terminal through a thin-film resistor, is 
used to control the circuit impedance at 3 GHz (the subharmonic of 
the 6-GHz band) without adding significant loss or mismatch at 6 
GHz.’ This was found to be necessary to eliminate frequency jumps 
during tuning, which occur when the subharmonic impedance is too 
high. 

At the end of the thin-film pattern (coaxial oscillator terminal), 
connection is made, using a bellows, to the center conductor of the 
coaxial line through the top half of the aluminum housing. 

The amplified rf signal reflected from the rwpatTr diode down the 
coaxial line is coupled by the center circulator to the output circulator. 
A small portion of the amplified output is capacitively coupled to the 
output detector circuit to provide the direct current for bay panel 
metering and alarm functions. This nondirectional coupling is approxi- 
mately 28 dB. The amplified signal passes through the output circu- 
lator, used as an isolator, and is coupled into a reduced-height wave- 
guide. Within the reduced-height waveguide, a waffle-iron filter® 
having a low-pass characteristic strips the amplified signal of any 
residual harmonic energy either generated by the IMpaTT or contained 
in the input signal. Following the waffle-iron filter, a four-step transi- 
tion couples the reduced-height waveguide to standard-height WR-159 
waveguide. 


3.2 Coaxial circuit 


A cross section of the coaxial line is shown in Fig. 7. At the bottom, 
just above the bellows contact to the thin-film circuit, is located a 
three-resonator, radial-line, band-stop filter? that is tuned to the 
12-GHz second harmonic of the 6-GHz common-carrier band. The 
filter prevents the second harmonic energy generated by the ImPATT 
diode from causing anomalous monitor circuit operation. Appropriate 
steps in the coaxial center conductor in the filter section provide a 
good match across the 6-GHz band. The center conductor tip is 
spring-loaded against the rmpatr diode, which is held centered at the 
upper end of the coaxial section. A large, finned heat sink contacting 
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Fig. 7—Cross section of the coaxial line section. 


the back of the diode provides diode cooling and eliminates the need 
for forced convection. 

IMPATT tuning is accomplished by a movable quarter-wavelength 
coaxial transformer and four capacitive power-adjustment screws 
located radially around the coaxial line at a point nominally an eighth- 
wavelength from the end of the transformer. The position of the trans- 
former relative to the impatTtT diode primarily determines the frequency 
of operation. The transformer is moved using a large-diameter knurled 
ring, shown just below the heat sink in Fig. 4. This ring is mechanically 
coupled to the transformer through two slots in the coaxial line. These 
slots are completely covered by the transformer and ring to prevent 
rf leakage. 

The transformer section characteristic impedance is designed to 
produce 1 watt at the amplifier’s output port (nearly 1 dB more at the 
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diode’s location) with the highest Q diode and with the power adjust- 
ment screws adjusted flush with the inner diameter of the coaxial 
outer conductor. For all other diodes that would give greater than 1 
watt with this transformer, capacitance is added using the four power- 
adjustment tuning screws. When transformed through the coaxial 
circuit to the diode position, the added capacitance appears as an 
increase in the resistive part of the circuit impedance. Increased circuit 
resistance reduces the generated power down to the 1-watt level, where 
optimum system performance occurs. This power adjustment is made 
on diodes having low values of Q; the increase in circuit Q because of 
the screw insertion is counterbalanced by the lower diode Q, so that 
the overall external Q is not increased by this power adjustment when 
compared with high Q diodes requiring little screw penetration. 


3.3 Circuit tuning 


The circuitry within the die-cast housing is initially tuned in the 
factory with a 7-mm precision connector located in place of the IMPATT 
diode and heat sink. During the initial tuning, the transformer is not 
installed and the power-adjustment screws are adjusted flush with the 
coaxial-line inner surface. Tuning of all ports of the three circulators 
to better than 30-dB return loss is accomplished across the 8-percent 
common-carrier band. The three ports of the center circulator are 
tuned over a slightly wider band to include the extremities of the 
locking bandwidth of amplifiers operated on the end channels of the 
common-carrier band. 

By matching the diode port of the center circulator to achieve this _ 
broadband high return loss, the oscillator circuit Q is essentially 
determined by that of the quarter-wavelength transformer and the 
short section of 50-ohm line from the transformer to the diode. In 
practice, Q..2 of the oscillator is determined largely by the mparr diode 
wafer. The transformer position and power adjustment screws permit 
adjustment in the field of any amplifier to 1 watt on any channel as- 
signment with any diode. The tmpaTt diode is replaceable in the field by 
simply removing the heat sink and inserting a new diode in the coaxial 
line against the spring-loaded center conductor. 

The amplifier cost has been kept low by the use of thin-film inte- 
gration, casting technology, and laser cutting of ceramic substrates. 


IV. AMPLIFIER PERFORMANCE 


Ten models were constructed in the laboratory, and information 
was conveyed to the Western Electric Company, who is now produc- 
ing the unit. Measurements of intermodulation distortion indicate 
that the distortion products are small and that system performance 
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can be accurately predicted on the basis of power output and FM 
thermal noise with no correction for distortion. Performance shown 
in Fig. 1 can be readily obtained using the silicon impart diodes in 
manufacture. A detailed evaluation has been completed of a Tm-2 
system in Ohio that includes eight factory-built impaTT amplifiers. 
Satisfactory operation was noted over a 10-month test period. 


Vv. SUMMARY 


A 1-watt, 6-GHz silicon impattT diode amplifier has been developed 
and is being manufactured for use as the transmitter power amplifier 
in short-haul radio systems. The amplifier operates with a nominal 
gain of 20 dB and a noise figure of less than 52 dB. The noise contribu- 
tion of the impaTtT amplifier is substantially thermal noise, with inter- 
modulation distortion negligible. The de-to-rf efficiency is 4 percent. 
The amplifier includes integrated input and output rf power monitors 
and harmonic suppression circuitry. 

The input monitor circuit furnishes the input information for the 
power-supply squelch circuit. If the input rf level drops low enough 
so that the locking bandwidth becomes small, the power supply is 
turned off, preventing the oscillator from free-running out of the 
assigned frequency range. The de power is automatically restored when 
the input level returns to normal. 

The low cost and reliability of this impatTr amplifier make it an 
attractive rf output device in short-haul applications. 
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Line-of-Sight Paths Over Random Terrain 


By E. N. GILBERT 
(Manuscript received September 5, 1974) 


Line-of-sight paths are important as VHF radio channels. In a mobile 
radio system, for example, the landscape determines the communication 
- possibilities in a complicated way. This paper analyzes a simple model of 
rough terrain to relate statistecal terrain properties to line-of-sight paths. 
The model is constructed from conical hills, all the same height, distributed 
at random over the surface of a spherical earth. 

The parameters of the. model are the earth’s radius a, the density o 
of hills, and the grade g of the hills. Although a simpler planar model is 
obtained by letting a—~, a finite spherical earth is needed for most 
questions. Assuming that a base station is located at the peak of a hill, 
the most interesting line-of-sight paths are those from a typical hilltop. 
A large number of statistics of these paths are then derived, usually as 
simple functions of a, c, and g. These include properties of paths to other 
peaks, to the horizon, and to random points on the ground. 


l. INTRODUCTION 


Very-high frequency radio propagation is often said to resemble 
optical propagation. A line-of-sight path provides a good radio channel; 
a path blocked by the terrain does not. With the aid of a topographic 
map, one can determine whether a path Q:Qz2 is a line-of-sight path. 
Essentially, one must plot the ground elevation profile along the path 
to see whether the ground intersects the straight line segment QiQ>. 
This calculation must include the effect of the earth’s curvature. 
Atmospheric refraction is also accounted for by changing the earth’s 
radius to a fictitious value. 

Having done the calculation for one path Q:, Qe, we learn little about 
other paths. The region covered by a transmitter at Qu, i.e., the set of 
points Q visible from Q:, would be found by plotting ground elevation 
profiles along views from Q: at every possible azimuth angle. This 
region might represent the coverage of a Tv station or of a base station 
in a mobile telephone system. 

This paper analyzes a statistical model to give insight into the way 
coverage regions depend on properties of the terrain. The parameters 
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of the model are the radius a of the earth, a density o of mountains 
(or hills) per unit area, and a grade (slope) g of these mountains. 
Many statistical properties of terrain and paths are then derived as 
functions of a, o, g. These properties are means, or in some cases 
distribution functions, of the random variables that appear in the 
INDEX. Line-of-sight paths from a typical mountain peak receive special 
attention because a peak is the most likely site for a base station. 
Although the exact formulas contain integrals with unwieldy trigono- 
metric integrands, most of these formulas may be replaced by simple 
expressions, to a very good approximation. The expected area visible 
from a peak and the expected number of peaks visible from a random 
point on the ground are more complicated quantities, leading to 
integrals that are evaluated numerically. 


INDEX 
Altitude—egs. (6), (7), (8), Table I, Fig. 6. 
Area blocking—egs. (12), (17) to (20), (23), Table III. 
Visible—eq. (43), Table VII. 
Within horizon—eq. (35). 
Number of peaks visible: 
From a peak—eq. (26), Table V. 
From a point on the ground—eq. (87). 
Range from a peak: 
To furthest visible peak—eq. (31), Table V. 
To horizon—egqs. (33), (34), Table VI. 
To random visible peak—egs. (25) to (29), Table IV. 
To random visible point on ground—egs. (89) to (42), 
Vig. 15. 
Slope—Table IT. 


The earth’s radius a is an important parameter of the model. 
Although a simpler planar model is obtained by letting a—~, the 
planar model is inadequate for most statistics of interest. 

With a and ¢@ fixed, the terrain becomes rougher as g increases. As a 
rule, the model predicts more long line-of-sight paths and larger 
expected visible area for rougher terrains. However, in mobile radio 
these long paths are more important as sources of interference than 
as useful channels. 


Il. THE MODEL 


The terrain model will use conical mountains distributed at random 
in a Poisson pattern over the surface of the earth. Begin with a sphere 
of radius a miles (a may be the true radius of the earth, or something 
larger if atmospheric refraction effects are to be taken into account). 
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Place points at random on the surface S of this sphere using a Poisson 
process with density o points per square mile. Each Poisson point will 
represent a mountain peak, and so the sphere of radius a will be called 
the peak sphere. 

Each Poisson point P will be associated with a mountain-shaped 
subset M(P) of the interior of the peak sphere. The subsets 1 (P), 
M(P:2), :-: for the various peaks will overlap. Take the union of all 
the subsets M(P) to represent the earth. 

The simplest shape for M(P) is the cone consisting of all rays from 
P making angle <@ with the inward-pointing normal to S. This cone 
has to be truncated to keep it from extending beyond the peak sphere 
in the direction antipodal to P. The surface of the cone is tangent to 
an inner sphere, concentric with the peak sphere and having radius 
asin 6. Take M(P) to be the inner sphere plus the part of the cone 
that lies between P and the inner sphere. Figure 1 shows M/(P) 
shaded. 

With this construction, the terrain consists of conical mountains, 
all having the same height and the same grade g = cot 6. There may 
also be flat places where the earth’s surface coincides with the inner 
sphere. A flat spot occurs at any point that lies further than ($7 — @) 
radians away from all Poisson points. Flat spots are rare, except when 
the parameters o, @ are chosen to produce widely separated mountains 
having very gentle slope. 
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Fig. 1—Construction of M(P). 


LINE-OF-SIGHT PATHS 737 


Figure 2 is an elevation contour map for a typical random terrain. 
Some unrealistic features of the model are evident. The conically 
shaped mountains have circular contour lines. The peaks are distrib- 
uted chaotically instead of being arranged in rows (mountain ranges). 
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Fig. 2—Contour map. The symbols - 8 + 6 — 4: 2 denote altitude levels ordered 
from the peak sphere downward. 
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Figure 3 is a plane cross section through the earth in the same model. 
This figure is more interesting for the present problem because the 
existence of a clear line-of-sight path between two points depends only 
on the shape of such a profile. Note that the elevation curve in Fig. 3 is 
composed of convex ares (hyperbolas) that join in the valleys between 
mountains. But, at least, the maxima in Fig. 3 have different heights. 
Figure 4 shows random terrain as seen from one of the peaks looking 
out toward the horizon. The nearest and furthest peaks shown have 
ranges of 6 and 150 miles. The parameters were picked to match a 
particular portion of the Alps for which a panoramic photograph was 
also available. The deficiencies of the model are less evident in this 
figure. The curvature of the earth makes it less obvious that all peaks 
have the same height. 

In real terrain, it is sometimes possible to see part of a mountain 
even though the mountain’s peak is obscured from view. That cannot 
happen in this model, as will now be proved. Suppose that the view 
of a peak P; is blocked when the eye is at #. Then the line segment 
P,E contains a blocking point B, belonging to another mountain 
M (Pz). Now consider any other point P of M(P1). P must lie on some 
line segment P17, where J belongs to the inner sphere. Figure 5 shows 
the triangle HP,J. The segments EP and Bel cross at some point B 
in the triangle. B belongs to the convex set M(P2) because Bz and I 
belong. Then B is a point of M(P2) blocking the view of P. 

By making a >, one obtains a planar model of random topography. 
The peak sphere S becomes a peak plane. At a point Q, the land surface 
lies below the peak plane a distance 


y = g Min ||P: ~ Ql (1) 


where the minimization is over all Poisson points P;. Replacing S by a 
plane simplifies the analysis considerably but, unfortunately, it 
produces a much less realistic model. If Fig. 4 has been drawn for a 
planar model, every peak P; would have been visible. Even worse, 
Section VIII shows that the expected area visible from a peak would 
be infinite. For that reason, the extra complication of a spherical earth 
is really necessary for some questions about line-of-sight paths. 


lil. PARAMETER ESTIMATION 


The two parameters o, g = cot 6 can be chosen to fit the model to 
terrain measurements. One might estimate the density o by counting 
peaks. A difficulty is that one must then decide how big a hill must be 
to be counted. Surely every bump on the landscape ought not to count 
as a peak. This decision is avoided by using statistical properties of 
the point Q lying below a random point g on the peak sphere. The 
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Fig. 4—View of horizon from a peak. 
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Fig. 5—M (P2) blocks all of M(P:) from view if it blocks Py. 


altitude and slope of the terrain at Q are two useful random variables. 
Both depend on the angle y = < POQ from Q to the nearest peak P 
(see Fig. 1). Since a circular cap of angle yo on the peak sphere has 
area 27a?(1 — cos yo), the distribution function for <POQ can be 
written immediately, 


Prob (y S yo) = 1 — exp [— 2za’o(1 — cos yo) ]. (2) 


There is no natural sea level in the model, and so it will be con- 
venient to specify the altitude at Q by giving the depth y, measured 
from S down to the land. If y = 4a — 6, then ground level coincides 
with the inner sphere, i.e., . 


y=al—siné), y2ar- 8. (3) 
For smaller angles y, 
y = all — sin 6/sin (y + 4) ], Os 7 < ir —- 48, (4) 


as is clear from Fig. 1. These formulas, together with the distribution 
(2) for y, determine the depth distribution, 


Prob {y S a[1 — sin 6/sin (y + 6) ]} 
= 1 — exp [—2ra’o(1 — cos y)], O0OSy<3r-—8 (5) 
Prob {y S$ a(i — sin é)} = 1. 


Although one can easily tabulate the distribution function for y by 
substituting numerical values of y into (5), the distribution function 
is easier to visualize in a limiting case. Since a is a large radius, let 
a—o in (5). As one might expect, the formulas tend toward the 
depth distribution function in the planar model, 


Prob {y S Y} = 1 — exp [—or(Y/g)?]. (6) 


In this limit, s and g enter the distribution only via a single length 
parameter o—?g, which is an index of altitude variability. Thus, altitude 
distribution data alone cannot be expected to supply good estimates of 
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both g and c. Some simpler statistics are the median, 


Median (y) = (a7 log, 29?/c)? 
= 0.469707, 


and the moments, 
E(y*) = TL + 3k) (g?/m0)*?. 


Particularly, the mean is 


g = Ely) = 30749 (7) 
and the standard deviation is 
[Var (y) ]? = [(a7 — 2-%)g?/o]? = 0.26830-%g. (8) 


It is also possible to obtain (6) as an exact result for a spherical 
model in which the shape of the mountains is only approximately 
conical. That entails a new choice of the set 1/(P) in Fig. 1. Define the 
new shape so that the depth becomes 


y = 2ga sin 37, (9) 


where again y is the angle to the peak. At P, M(P) comes to a point 
approximating a cone of slope g. At the antipode to P, M(P) has depth 
2ga; then this model requires g < 3. Now the depth distribution for all 
y is again given by (5) but with the left-hand side replaced by 
Prob {y S 2ga sin $y}. But that is (6), exactly. 

For many values of g and o, the planar approximation (6) to the 
depth distribution (5) is very good. For example, Table I compares 
the planar approximation with some distributions having a = 3959 mi, 
the earth’s radius. In the table, the cones have grades g = 0.05, 0.1, 
and 0.2 and the density o is adjusted to fix the standard deviation in 
(8) at 528 ft (0.1 mi). Table I gives percentiles of the distribution as 


Table |— Altitude percentiles (in feet) 





Spherical Model Planar 
¢= 0.05 0.1 0.2 2 
a 0.0171 0.0683 0.2732 | 7/9 = 6.831 
0.1% ~1899 ~1964 1980 —1986 
1% _ 1378 ~1421 ~ 1439 ~ 1436 
10% —~691 —712 ~717 _719 
25% —315 327 331 — 339 
50% 70 63 62 61 
75% 402 400 399 399 
90% 642 641 640 640 
99% 896 896 896 896 
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altitudes measured upward from a common level, corresponding to the 
depth g in (7). 

Figure 6 is an altitude distribution for northern New Jersey. It 
was obtained from a topographic map by reading altitudes at 52 
points, 10 km apart in a rectangular grid covering latitudes 40°30’ 
to 41° and longitudes west of 74°. The altitudes ranged from 0 to 
1100 ft. Data for parts of New Jersey further south were not used; the 
topography of New Jersey is too variable for both north and south 
to be well represented by a single simple model. The planar model 
fits the observed points well, except at low altitudes. As an alternative, 
use the spherical model with a = 3959 mi. By taking g = 0.011, one 
obtains a maximum depth (3) near 1100 ft, so that low altitudes can 
be regarded as occurring on the inner sphere. Then o remains as a 
parameter to adjust for a good fit. 

The parameter g = cot @ is the grade at mountain peaks. At the 
random point Q, at angle y away from a peak in Fig. 1, the grade is 
smaller because the normal to the conical surface makes an angle 
iq — 0 — y with the vertical direction OQ. Thus, the grade at Q is 


,__ joot@+y) = (g—tany)/i+tgtany), if @6+y Sam 
9 0 otherwise. 


PROBABILITY 





0 100 200 300 400 500 600 700 800 900 1000 1100 
ALTITUDE IN FEET ABOVE SEA LEVEL 


Fig. 6—Altitude distribution for northern New Jersey. oe is for planar model 
with peak sphere at 1130-ft altitude and [Var (y) ]! = 400 
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Table Il — Percentiles of G/g 


g = small 0.1 0.2 0.5 
0% — 1.0000 — 1.0000 — 1.0000 — 1.0000 
1% —0.9995 —0.9995 —0.9995 —0.9994 

10% —0.9510 — 0.9506 —0.9492 — 0.9399 

25% — 0.7070 — 0.7053 —0.7001 — 0.6667 

50% 0 0 0 0 

75% 0.7070 0.7053 0.7001 0.6667 

90% 0.9510 0.9506 0.9492 0.9399 

99% 0.9995 0.9995 0.9995 0.9994 

100% 1.0000 1.0000 1.0000 1.0000 





The grade 0 occurs on the inner sphere. This result, together with (2), 
determines the distribution of the grade g’. In most cases, the grade g’ 
has high probability of being close to g; one should not expect this 
distribution to fit observed grade data well. 

At g, one might move in a random direction and ask for the slope G 
along the random path through Q. The slope, which depends on the 
angle ¢ between the path direction and the uphill direction, lies in the 
range —g’ SG Sg’. With some simple geometry, one finds 


G = g' cos g/(1 + g” sin? ¢)}. 


By using the known distribution for g’ and assuming a flat distribution 
for y, one can obtain a distribution function for G, This would be the 
distribution of the slopes G seen in cross sections like Fig. 3. A simple 
distribution is obtained only in the planar model limit, for which 
g =g= cot 6 identically : 


Prob {G S tan X} = 1 — w' arc cos {sin X/cos 6}. 


Table II gives the slope distribution in the planar model for several 
values of g. In the limit of small 9; the distribution function for G/g 
tends to 1 — a arc cos G/9. 


IV. BLOCKING REGIONS 


Suppose two points Q1, Qe are given, representing the positions of 
two antennas. In general, Qi, Qe can lie anywhere above the inner 
sphere. A clear line-o f-sight path exists between Q1 and Q2 as long as 
the straight-line segment Q1Q2 does not intersect any of the sets 
M (P;). The blocking region for Q1, Qo is the (open) set of points P on S 
such that Q:Q2 intersects M(P). The area of the blocking region enters 
into the probability that a line-of-sight path QiQ2 exists. The advantage 
of the conical mountains M(P) is that blocking regions assume simple 
shapes. 
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The simplest blocking region is one for a pair of points Q1, Q2 both 
on S. If the line Q:Q2 intersects the inner sphere, all 1f(P;) block the 
path. The blocking region then consists of the entire sphere S. If 
Q1Q2 misses the inner sphere, then blocking occurs at a point Q on the 
path Q:Q.2 if a peak P; lies too close to Q. If the depth of Q is y, then 
(4) gives the angle y to peaks P such that Q lies on the surface of 
M(P). Then blocking occurs at Q if a circular cap of angular radius 
y contains a peak. The pole of this cap is the radial projection q of Q 
onto S. The blocking region for the path QiQ. is the union of all the 
blocking caps for points Q on the path. These caps are largest midway 
between Q; and Qo, shrinking to points at Q: and Qs. Then the blocking 
region is lens-shaped, as in Fig. 7. 

Figure 7 shows two ares K, K’ which form the boundary of the 
blocking region. The argument that follows shows that K, K’ are 
actually ares of circles. Figure 8 is another view of the peak sphere 
projected directly along the line Qi, Qo. Two planes, 7 and x’, can be 
drawn through Q:, Q2 and tangent to the inner sphere, say at C and C’. 
These planes project to lines in Fig. 8. The planes 7 and 7’ intersect S 
in two circles, centered at C and C’ and both passing through Q: and 
Q:. Since M(P) is the convex hull of P and the inner sphere, 7 is a 
supporting plane of M(P) as long as P lies below 7m (i.e., in the half- 
space containing the inner sphere). Then M(P) does not block the 
path Q:Q.2 if P lies below z, or below zw’. The part of S lying above 
both 7 and z’ appears shaded in Fig. 8. Suppose P belongs to the shaded 
region. Project the triangle C’PC, a subset of /(P), onto the plane 
of Fig. 8. The path Q1:Q2 projects to a point lying inside this projected 
triangle. Then Q:Q2, a chord of S, must intersect the triangle C’PC. 





Fig. 7—Blocking region for two points Qi, Q2 on the peak sphere. 


746 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1975 


INNER 
SPHERE 





Fig. 8—Another view of the blocking region. 


The point of intersection is a point of M/(P), which blocks the path. 
Thus, the shaded region, bounded by ares of the circles S () 7, Sf) 7’; 
is the blocking region for Q1Q2. 

The area A(Q:, Q2) of the blocking region in Fig. 7 will now be 
expressed as a function of the angle 2o = 2Q.:0Q». Project the centers 
C, C’ in Fig. 8 radially out to c, c’ on S. Figure 9 is another view of S 
showing c, c’ as the poles of two circular caps bounded by S () 7 and 
S) 7’. The angular radius of both caps is $7 — 6, as is clear from 
Fig. 8. The chord Q:Q2 subtends some angle 2a = /Q:cQ2 at c. Using 
the spherical sine law in the right triangle Qi, c, $(Q1 + Q2), one may 
determine a from 

sin a = sin p/cos #. (10) 


The cap with pole c has area 27a?(1 — sin 6) and the sector included 
within angle 2a has area 


A, = 2aa?(1 — sin 6). 
Also, the triangle QicQ2 has area 
Ar= (2a + 26 — 1) a’, 


where 8 = Q1Q:C1 = Q2Qic. The sine law may be applied to triangle 
QicQe to find B 
sin 8 = cos 6 sin 2a/sin 2p. (11) 
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Fig. 9—Angles used in deriving area A of blocking region. 


The difference As — Ar is 4A (Qi, Q2). Thus, 
A(Q:, Qo) = a?(2r — 48 — 4a sin 6), (12) 


where (10), (11) give a, B. 

The blocking region is more complicated if Q1, Qe or both are not 
on S. As in Fig. 7, each point Q on Q:Q, is blocked by peaks lying in a 
circular cap of radius y given by (4); the blocking region is the union 
of these caps. Let Q:, Q2 be the points where the extended line Q:Q> 
meets S. The blocking region for Q1Q> is a subset of the blocking 
region for QiQ2. As shown in Fig. 10, the blocking region consists of 
the caps for blocking at Qi and Q: plus the part of the blocking region 
for QQ. that lies between these caps. The centers of the two end caps 
are the points q1, ge obtained by projecting Q:1, Q2 radially onto S. 

The two end caps have a special role in the blocking. Normally, 
Q1, Q2 are known to lie above ground, and so the two end caps are 
known to contain no peaks. If the ground levels below qi, q2 are known, 
then peaks Pi, P2 must exist somewhere at the appropriate angles 71, 
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yo away from qi, ge. In Fig. 10, Qi, Q2 are assumed to be at ground 
level; then Pi, P2 lie on the boundaries of the end caps. The mountains 
M(P:), M(P2) on which Q:, Q2 lie can themselves block the path QQ. 
Thus, in Fig. 10, P2 blocks the path because it lies in the blocking 
region. To compute the conditional blocking probability for the 
configuration in Fig. 10, one must know both the area of the part of 
the blocking region that lies outside the end caps and also angles ¢1, ¢2 
that limit where Pi, P2 can lie to cause blocking. In the applications 
that follow, it will suffice to let Q: lie at a peak Qi = P; and let Q2 bea 
point at ground level. That simplifies Fig. 10 to Fig. 11. 

Let 22 = ZP,10q2, the angular distance along the are Pigq2. The depth 
at Q» determines the angle y2 of the end cap. The sides of the spherical 
triangle Pigec are now known, and so its anglesB = /@q2Pic, 
C2 = ZqcP1, tr — ge = ZPigec are determined. One finds 


1 — cos (22 — ye) + 9 sin (22 — y2) |? 





tan$g: = | 

onan g Sin (22 + y2) — 1 + cos (22 + Y2) o 
sin 8 = sin g2 cos (6 + y2)/cos 6 (14) 
sin f2 = Sin ¢g2 Sin Z2/cos 0. (15) 


The blocking area is twice the area of the half of the blocking region 
above the line PiQ2 in Fig. 11. That half can be obtained as a sum of 
two parts. One part is a sector of angle 7 — ge from the end cap; its 
area is (t — g2)(1 — cos y2)a?. The other part is obtained by removing 


c! 
° 





Fig. 10—Blocking region for points Q1, Q2 which are not peaks. 


LINE-OF-SIGHT PATHS 749 


Q,=P, + 








fe) 
Fig. 11—Blocking region for Q: = P, a peak, Q2 not a peak. 


the triangle Piqoc of area (8 + £2 — ge)a? from a sector centered at c. 
The sector has area ¢2(1 — sin @)a?. These areas may be combined to 
express the area of the blocking region in the form 


A(Q1, Qe) = Ao i Ag, 


where - 
Ag = 2r(1 — cos y2)a? (16) 


and 
Ao = (2¢2 cos y2 — 28 — 2f2 sin 6)a?. (17) 


Ag is the area of the end cap and Ais the area of the remainder of the 
blocking region. 
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Although the blocking area A(Qi, Q2) for the general situation of 
Fig. 9 will not be needed, it can be obtained in the form 


A(Q1, Q2) = A(Q1, Q2) + A(Qi, Q2) — A (Qi, Q2). (18) 


Note that formulas like (16) and (17) give A(Q1, Q2), 4(Q:, Q2) while 
(10), (11), and (12) give A(Q:, Q2). Likewise, with a change of sub- 
scripts, (13) gives ¢1 as well as go. 

These formulas can now be used to obtain the path probability 
p(Q1, Q2), the conditional probability that a clear line-of-sight path 
exists between given points Qi, Qe. When Qu, Qz are on S, as in Fig. 7, 
p(Q1, Qe) is just the probability that the shaded region of area A (Q1, Q2) 
contains no peaks. Then 


p(Q,, Qe) = exp [—oA (Qi, Qe) ], (19) 


with A(Q:, Qe) given by (12). 

The situation in Fig. 11 is more complicated. Q: = Pi, a peak, and 
Qs is supposed to lie on the ground. Then the cap of area Az is known 
to be empty. Two conditions must hold if the entire blocking region 
is to be empty. One is that the peak P2 of the mountain on which Q2 
lies causes no blocking. Since P2 is equally likely to be anywhere on 
the boundary of the cap around Qs, there is probability 1 — g2/a that 
P, does not block the path. The second condition is that the remainder 
of area Ag of the blocking region is empty. Then 


P(Q1, Q2) = (1 — ¢2/m) exp (—oAd), (20) 


where (13) and (17) give g and Ao. Formula (20) applies as long as 
v2 < 2. It is also possible to have yz = 22. In that case, Q2 lies on the 
mountain 1/(P1); the path in question runs from the peak Qi = Pi 
to Qe along the surface of the cone M(P;1). Whether or not such a path 
is to be considered blocked is a matter of definition. Here M(P1:) is 
regarded as an open set so that the path is not blocked. As y2 — 2s, 
one finds g.—>0 and Ay—0 so that p(Q:, Q2) > 1, ie., (20) con- 
tinues to give the correct probability in the limit. 

Another limiting situation, y2 — 0, illustrates an important distinc- 
tion between Figs. 11 and 7. In the limit Az — 0 and Ay becomes the 
area of a lens-shaped region, such as shown in Fig. 7. Then the ex- 
ponential factor in (20) becomes the path probability (19) for the 
two peaks Qi(=P;), and lim Q2. However, (20) contains an extra 
factor (1 — ¢g2/r) which approaches 3, not 1. This disagreement 
between (19) and (20) is explained as follows. From a point Q2 near a 
peak P:, one can look over a 180-degree view; the mountain M(P2) 
blocks the other 180 degrees. Then the factor 4 in (20) is needed to 
account for possible blocking by M(P:2). But exactly at the peak 
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Q. = Ps, the mountain 1/(P2) no longer interferes in any direction. 
Then no factor 3 is needed in (19). The discontinuity in p(Q1, Q2) as 
Q2 — Pz could be avoided by assuming that the antenna location Q» 
lies at some known positive height h above ground. 


V. PATHS BETWEEN PEAKS 


The simplest blocking regions were for paths QiQ2 with both end 
points on mountain peaks. The path probability p(Q1, Qe) in (19) can 
now be used to derive some interesting properties of peak-to-peak 
paths. In this section, Q: will be a given peak Py. Qe will be another 
peak selected at random. An element of area dA(Q2) on the peak 
sphere S has probability odA(Q2) of containing a peak Qs. Then 
op(Q1, Q2)dA (Q2) is the probability that the element contains a peak 
Qe which is visible from Q,. 

Let d(Q:, Q2) denote great circle distance between Q: and Qe. Let 
~,(d) denote the random variable which is the sum 


2k (d) = x d(Q:, P;)* (21) 


of kth powers of distances from Q, to all other visible peaks P; lying 
within distance d[d(Q:, P:) S d]. The element dA(Q2) contributes a 
term d(Q:, Q2)* to 2.(d) with probability cp(Q:, Q2)dA (Qz2). Thus, the 
expected value of 2;(d) is 


HEI] = [ f dQ, Q)*p(Qs, Q)4A(Q), (22) 


where the integral extends over all points Q2 in the cap d(Q1, Q2) S d. 
Another random variable 2, is a sum like (21) extended over all 
visible peaks, at any distance from Q;. The mean H(2;) is an integral 
(22) over the entire sphere. Evaluating (22) will give the mean number 


Table il] Blocking area A(Q,, Q2) in square miles, as given by 
exact and approximate formulas (12) and (23) with 
range d(Q,, Qz) = 100 miles 


Grade g Exact Approximate 
0.01263 7887.7 3333.1 
0.02 2430.9 2104.9 
0.05 857.6 842.0 
0.1 423.2 421.0 
0.2 210.2 210.5 
0.5 84.1 $4.2 
1.0 42.0 42.1 
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of visible peaks E(2 o) and other information about the distances to 
visible peaks. That could be done numerically, using (10), (11), (12), 
and (19). However, the approximation that follows simplifies the 
evaluation. 

The approximation is one which holds when a is so large that angles 
2p between visible peaks can be considered small. The planar model 
has A(Q1, Q2) = 0 and p(Q:, Q2) = 1, which is too rough to make 
sense in (22). Instead, the first nonzero term in a series for A (Qi, Q2) 
in powers of p will be used. Expansion of the exact formulas (10), (11), 
and (12) is laborious but straightforward: 


a = p(1 + g*)?/g + p®(1 + g?)3/(69%) + 0(p°) 
B = 3m — p/g — (2g? + 1)p?/(6g3) + 0(p°) 
A(Q1, Qe) = a?(2p)?/(6g) + 0(p'). 
= r3/(6ga) + -°-, 


where r = 2ap is the great circle distance d(Qi, Qe). 

For a simpler, more intuitive, derivation, one may find the size of 
the circle about a typical point gq along the path Q:, Qe in Fig. 7. If z is 
the great circle distance from Q: to q, then the ground level below q 
lies at depth y satisfying 


(a — y) cos (2 — 4r)/a = acosar/a 
or 


y = 2(r — 2)/(2a) + +, 


The radius ay of the circle about q is approximately y/g, and the 
blocking area is approximately 


A(Q1,Q) = [" ade 


a a(r — z)dz/(2ag) 
r/(6ga) + ---, (23) 


I 


A (Q:, Q2) 


as before. 

From the form of the series used in deriving (23), one may predict 
that the rate of convergence is determined by the ratio p/g. Table III 
shows that (23) does give a better approximation for large g than for 
small. In Table III, a = 3959 mi; a large range, 100 mi, was used for 
a severe test of the approximation. At grades g smaller than 0.01263, 
a 100-mi path between peaks is blocked by the inner sphere. The small 
p/g condition is another way of requiring that the path Q1Q:2 clears 
the inner sphere by a safe margin. 
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Now use the approximation (23) for A(Q1, Q2) in (19) to evaluate 
the integral (22) for the expectation [2;(d)]: 


d/ (2a) 
E(2;.(d)] = 2x0 i (2ap)**! exp { —4a?cp?/3g} 2adp 
0 


(d/D)4 
E[z:(d)] = (240/3)D*? I u®D2 exp (—u)du, (24) 
0 
where 
D = (6ag/c)}. 


The integral in (24) may be expressed in terms of the incomplete 
gamma function, 


E(2:(d)] = (2010/3) D'? {TL (k + 2)/3] — E(k + 2)/3], (d/D)%}, 
or the X? distribution function, 
E[2i(d)] = (2940/3) D* YT (k + 2)/3]P(x?|r), 
where 
x2 = 2(d/D)8 
and the number of degrees of freedom is 
= 2(k + 2)/8. 


Although the approximation (23) becomes poor at long ranges, the 
integrand is very small there. Thus, (24) can be expected to hold even 
for long ranges. In particular, the expectation #(2;), for visible peaks 
of all ranges, may be approximated by letting d > in (24): 


E(2x) = (240 /3)D* PT (k + 2)/3 ]. (25) 
For the special value k = 0, (25) gives the mean number of peaks 


visible from Q:: 


E(2o) = roD*T (5/38) 
9.3645 (a2g2c)? 


= 2344(g’c)? if a = 3959 mi. (26) 


I 


Note, as predicted earlier, that the mean number of visible peaks 

tends to infinity in the limit of large a@ (planar model). When k = 1, 
(25) simplifies to 

E(21) = 4rag. (27) 

As o increases, (26) shows that the mean number of visible peaks 

increases, but (27) shows that the mean sum of distances to visible 


peaks remains unchanged. This indicates that visible peaks tend to 
be closer for large o than for small. One way to define a range for a 
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“typical”? peak is to form the ratio 


D, = E(21)/E(2o) = D/T (2/8) 
= 0.738487D 
= 1,34190(ag/o)!. (28) 


VI. RANGES BETWEEN PEAKS 


One might ask for a probability distribution for the range d from 
a peak P; to a randomly chosen visible peak P # P;. The random 
process for choosing a peak must be specified with care. Perhaps the 
most natural process would be this. Construct a random landscape 
and choose a peak P from the set of 2o visible peaks, all peaks equally 
likely. Then ask for the probability that P is one of the 2o(d) peaks 
within range d of P;. Given a landscape, the conditional probability 
that P is within range is 2o(d)/Zo. Then the unconditioned probability 
is E[Z(d)/Zo]. Unfortunately, the expectation is hard to obtain [there 
is also a question of giving an appropriate meaning to 2o(d)/2» when 

By using a different random process, one obtains a simpler distri- 
bution. Construct a trial random landscape and pick one of the peaks 
P at random, this time from the set of all peaks on the entire sphere S. 
P may not be visible. If not, discard that trial and construct a new 
landscape. Continue constructing landscapes and choosing peaks until 
the chosen point P is visible. Then ask for the probability p(d) that 
P lies within range d. 

To determine p(d), note that the total number of peaks on the entire 
sphere has the Poisson distribution with mean 47a?c. The argument to 
follow assumes that this number is large, so that the number of peaks 
actually obtained is almost always very close to its mean value. Then 
the probability that P is visible is H(2o)/(42a?c). If q[ Zo, Zo(d) | is the 
joint probability for 2» and Zo(d) in each trial, then g[Zo, 2o(d) ]Z0/ 
(47a?c) is the probability that a trial has 2» visible peaks, 2o(d) within 
range d, and that P is visible. The joint probability for the numbers Zo, 
Yo(d) of the landscape, selected when P is visible, is q’[2o, Zo(d) ] 
= q[2o, Zo(d) |Z0/H (Zo). The probability that P lies within range d 
is obtained as a sum over Yo(d) and Lo 


p(d) = X q'[2o, Zo(d) J20(d)/Zo 
E[2o(d) J/E (Zo) 

p(@) = 1 — T[2/8, (d/D)*]/T (2/8), (29) 
the last line following from (24). 


It is clear from this derivation that the second random process 
tends to select random landscapes with larger 2» than the landscapes 


I 
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Table IV — Probability p(d) = E(3o(d))/E(%o) that a randomly 
chosen visible peak lies within range d 


d/D Probability 
0.25 0.06880 
0.30211 0.1 
0.48595 0.25 
0.5 0.26361 
0.72212 0.5 
0.75 0.53050 
1 0.77518 
1.20507 0.9 
1.82182 0.95 
1.5 0.98440 
1.55886 0.99 
1.81350 0.999 
2 0.99983 


of the first process. However, 2» may be expected to have a highly 
peaked distribution, in which case 2» is nearly always close to E(2»). 
Then q(-, -) and q’(-, -) are nearly the same, and (29) is also a good 
approximation to the range distribution for the first random process. 

Equation (29) provides numerical values for the range distribution 
in Table IV. 

Another random variable of interest is the range to the furthest 
visible peak. Even the expectation of this maximum range is hard to 
derive. However, a simpler ‘typical’? maximum range is the range dm 
such that the expected number of visible peaks with ranges d > dm is 
just 4. Then d, satisfies 


E (Zo) — ELZo(dn)] = 2, (30) 
and (24) shows that 


I * actexp (—u)du = 3/(4neD?). (31) 
(dm/D)* 


Table V— Mean number of visible peaks E(%o) and range d,, 
such that E(3o-S0(dm)) = 


oD? aga E(2o) din /D 
3.11 0.84 8.8 1.30 
5.69 5.12 16.1 1.40 

11.3 40.1 32.1 1.50 

24.5 409 69.5 1.60 

58.4 5528 165.6 1.70 

153.6 100572 435.5 1.80 

450.5 2.54 X 10° 1278 1.90 
1477 8.95 X 107 4189 2.00 
5447 4.49 X 109 15448 2.10 
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Table V gives numerical values of d,/D as a function of oD? 
= (36a7g’c)*. H(2Z»), which also depends on oD? as shown by (26), 
also appears in the table. Note that d,, is not just a function of a single 
product of powers of a, g, «; it has a more complicated form (ag/c)? 
X function (a’g?c). 

The integral in (31) is a rapidly decreasing function of d,/D. Then 
the numbers in Table V would not change much if d, were redefined 
with the term 4 in (30) replaced by any other number of the same order 
of magnitude. For the same reason, d,, can be expected to be a good 
approximation to the mean range to the furthest peak. 


Vil. THE HORIZON 


The approximation (23) will now be used to derive properties of 
the range from a peak P; to the horizon at a random azimuth angle. 
The range to the horizon is a more interesting random variable than 
the range to a random visible peak. As has been noted, it is not always 
clear what to count as a peak in a real landscape. But the horizon 
has no ambiguity. 

Look from P with a fixed azimuth angle. One sees only sky at high 
elevations and ground at low elevations. The horizon point is the limit- 
ing point at ground level which has the highest elevation angle. The 
distance z from P; to the horizon is the range of interest here. 

Figure 11 will now be used to derive the conditions under which Q» 
is the horizon point, as seen from a peak P. If Qe is the horizon point, 
the entire straight line path P1Q; in Fig. 11 must intersect the ground 
only at Qe. Then the entire lens-shaped blocking region for Q2 must 
contain no peaks. But the depth at Q2 determines the circle on which 
a peak must lie. This circle appears in Fig. 11 inside the (open) blocking 
region. The only place that this peak can be now is on the boundary 
of the blocking region at one of the points of tangency 7, 7’. 

Figure 11 shows the usual situation in which the horizon point is 
not on the inner sphere. There is small probability that Q2 is on the 
inner sphere. In that case, go is at angle $7 — 6 away from P,, the 
centers c, c’ coincide with q2, and the blocking region is bounded by the 
circle through P, with center gz. There is no second peak on the bound- 
ary of the blocking region; Q2 lies on M(P). 

To find the probability distribution for the horizon range, one may 
first find the joint distribution for that range and the range r to the 
intersection point Q:. Since r is the largest range for which the corre- 
sponding blocking region is empty, the probability distribution func- 
tion for r is P(r) = 1 — exp [—oA (Pi, Q:) ]. In this derivation, the 
approximation (23) will be used to write 


P(r) = 1 — exp [—(r/D)*]. 
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Figure 12 shows Q2 lying at a range between r and r + dr, an event 
of probability dP(r). Given this position for Q2, the band between the 
boundaries of the blocking regions at r and r + dr contains the peak 
on which the horizon point lies. The shaded part of this band is the 
region in which the peak must lie so that the horizon point Q2 will 
have range 22 S z. The conditional probability function for the horizon 
range is just the ratio of the area of the shaded part of the band to the 
total area of the band. To simplify that calculation, one may replace 
the dotted line by a great circle that crosses P:P2 at right angles. That 
approximation leads ultimately to the conditional distribution 


Prob {horizon range S z|r} = (z2/r)?, OSzsSr. (32) 


The details are omitted because the result can almost be guessed 


immediately from the roughly triangular shapes of the two parts of the 
shaded region. 


Now the unconditional probability distribution for the horizon range 
is obtainable from (32) by integrating 


Prob {horizon range S z} 
= [apa + f° emar@) 
= (z/D)T[3, (z/D)*] + 1 — exp [—(z/D)*]. (83) 


Table VI gives numerical values. 


r+dr 





Fig. 12—Horizon at range S z. 
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Table VI — Distribution function for range z to the horizon 
looking from peak P, with a random azimuth angle 


2/D Probability 
0.21417 0.1 
0.25 0.13625 
0.35618 0.25 
0.5 0.42355 
0.56305 0.5 
0.75 0.70432 
0.79977 0.75 
1.0 0.88853 
1.02324 0.9 
1.15749 0.95 
1.25 0.97109 
1.40527 0.99 
1.5 0.995247 
1.67110 0.999 


The moments of the horizon range z are easy to find. From (83), the 
probability density for z is 


22 / * 2dP(r). 
The kth moment of z is 


E (2) 


2 i . ie / 7 r—dP (r)dz 
[2/(k + 2) ]E(r*). 


The last line is obtained by integrating by parts. The expectation on 
the right is another integral that can be evaluated in the manner of 
(24) and (25). The final result is 


E(2') = 2D'T(1 + (4/3) 1/(k + 2). (34) 


Equation (34) with k = 2 is particularly interesting. If z(¢) is the 
range to the horizon when looking with azimuth yg, then the mean 


area within horizon range is 
Qa 
a(3 [ A(¢)de) 
2 Jo 


= rE(z) 
= 4qI'(5/3)D? 
= 1,41803D°. (35) 


E(area within horizon) 


This expectation is only an upper bound on the mean area visible. 
For, as is clear in Fig. 4, there are points within horizon range that are 
obscured from view. 
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It is interesting to compare Tables IV and VI. At any given proba- 
bility level, the range to the horizon is smaller than the range to a 
randomly chosen visible peak. This may be surprising at first. How- 
ever, each visible peak is itself a point on the horizon. As seen in Fig. 4, 
the horizon consists entirely of small line segments extending down the 
sides of the mountains from the visible peaks. The line segments for 
distant peaks tend to subtend smaller azimuth angles than the seg- 
ments for nearby peaks. Picking an azimuth at random, one is more 
likely to find the horizon point on one of the nearby visible peaks 
than on one far away. 

Another expectation that exhibits the same effect is the mean 
azimuth angle between the horizon point and the peak of the mountain 
on which the horizon point lies. The ranges z and r determine this 
angle. Without belaboring the details, one can approximate this angle 
by its tangent and make the further approximations by which (83) 
was derived. The expected angle is found to be 


E(angle) = EL(r — z)/(2ag)] 
= T'(4/3)D/(6ag). 
That result can be stated in a more illuminating way: 
E(%o)E(angle) = 7I'(4/3)T (5/3) 
= (273-%) (27) 
= 0.40306 (27). 


By contrast, if E(2o) peaks were evenly distributed in azimuth with 
angular separation 27/H(2o), one would obtain 


E(2o)H (angle) = 0.25(27). 


The larger factor 0.40306 again occurs because of the variability of the 
angles which visible mountains subtend on the horizon. 


Vill. COVERAGE AREA 


The coverage set for a point P is the set of points Q such that a line-of- 
sight path PQ exists. In vHF radio applications, the coverage set of P 
is the set of points Q to which an antenna at P can radiate a strong 
signal. This section will estimate the mean coverage area C, the expected 
area of the coverage set for P = P1, a peak. 

The coverage set can have a very complicated shape. Figure 13 
shows one coverage set. In Fig. 13, the peaks are not in a Poisson 
pattern; to simplify the drawing, the peaks were located at points of a 
regular lattice. The coverage set contains the entire mountain on 
which P lies plus parts of adjacent mountains. These points alone 


760 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1975 








O 


Fig. 13—A coverage set. 


constitute a hexagon of area 4/c. In addition, the coverage set contains 
many smaller isolated patches on more remote mountains. These 
small patches can be so numerous that they represent most of the 
coverage area. 

If Fig. 13 represented the coverage set of a base station in a mobile 
radio telephone system, the station would only serve the hexagon of 
area 4/c. The other small patches would lie in the service areas of 
other stations, and so these patches would represent places where the 
given station can interfere with other stations. 

As in Fig. 11, let Pi be a given peak and Q» another point at ground 
level. Suppose the distance r from Pi to Q2 is known, i.e., the angle 
z=r/a in Fig. 11 is given. Let f(r) denote the probability that a 
line-of-sight path P1Q» exists. Since an element of area dA (Q2) at Q» 
belongs to the coverage set of P: with probability f(r)dA(Q2), the 
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mean area covered is 


Ge / | f(r)dA (Qs) = 2a? i " f(r) sin zedzo. (36) 


Before attempting to evaluate f(r) and C, the integral (36) will be 
given a second interpretation. Now consider Qe at a fixed location 
and count the number of peaks visible from Qe. The probability that 
an element of area dA(P:) contains a peak P, visible from Q: is 
of(r)dA(P). Then the mean number of peaks visible from Qz is 


E (visible peaks) = o / / f(r)dA(P) 


2rare / / f(r) sin zedz2 
= aC. (37) 


I 


Equation (37) can be used to derive very simple bounds on C. 
Clearly, more peaks are visible from a point Q2 at high altitudes than 
at low. If Qe were itself a peak, the mean number of visible peaks 
would be H(2o), given by (26). But Q2 has probability zero of being 
exactly at a peak. If Q2 is at any slightly lower altitude, Q2 is on the 
side of a hill which obscures 180 degrees of the view from Qe. Thus, 
E(visible peaks) S $E(2o), and (37) shows 


CS E(2o)/(2c). (38) 


Curiously, the right-hand side of (38) is exactly the mean area within 
the horizon as given by (35). Then (88) is a bound that was obtained 
in Section VII. 

At the other extreme, Q2 might be on the inner sphere, where no 
peak is visible. In most cases, that event will be so unlikely that it 
will be safe to say that the worst reasonable possibility is that Qe is 
down in a valley near the point where three mountains meet. Here, 
the three mountain peaks are visible and so one concludes C 2 3/c. 

To obtain f(r), and hence C, recall that (20) is a formula for the 
path probability p(Q:1, Q2), depending on the altitude y at Qe. To get 
f(r) one may average p(Q1, Q2) over y (or yz). This average may be 
expressed as a sum of two terms which account for the possibilities 
that Q2 belongs to the same mountain M(P;) as Q; or to a different 
mountain. 

In Fig. 11, if y2 = zz, then Qe lies on the mountain M(P,). This 
event has probability exp {—2mca?(1 — cos z2)}. The path probability 
p(Q1, Q@2) = lifr S (gm — Oa. lf r > ($m — 6)a, then Q: lies on the 
inner sphere, the path P1Q. is blocked, and p(Q:, Q2) = 0. 
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If yo < zo, then Q» lies on a different mountain M(P2). This possi- 
bility contributes a second term to f(r), 


J", P(@s, Q2)at1 — exp (—0s)}, 


where Az and p(Q1, Qe) are given by (16) and (20). 
For r S (4m — 6)a, the two terms combine into 


f(r) = exp {—22ca?(1 — cos 22)} 
a ae (1 — g2/m) exp l—o(Ao + Az) ]odAx, (89) 


where (16) and (17) provide Ae, Ao. A similar formula applies when r 
is larger, but f(r) is very small at such large ranges. 

One could find f(r) to any desired accuracy by evaluating the integral 
in (39) numerically. Instead, (39) will be replaced by a simpler approxi- 
mate formula. Since a is large, the first term of (39) is approximately 
exp (—o7r?); also, Az = wx? where x = ay2. The approximations to ¢ 
and Ao which follow are not uniformly good but are intended to apply 
in situations that contribute most of the coverage area C’. Except at 
very short ranges r, a typical blocking region is more elongated than 
that shown in Fig. 11. Figure 14 is more typical. Then g2 = $2, 
approximately. With that approximation, the shaded region in Fig. 14 
has area Ay + $A2. It consists of a triangle, of area xr, and two extra 
lens-shaped pieces. The two extra pieces can fit together into one lens 
of exactly the shape of the blocking region for two peaks at separation 
r. Then the two extra pieces combined have area r*/(6ga), as in (23). 
Now the exponent «(Ao + Az) in (39) is approximately (r/D)* + $2027. 

To substitute these approximations into (39), write 


Y = onr 


X = r°/(6ga) = (r/D)* = [1(5/3) Y/E (Zo) }} oy 





Fig. 14—Approximation of Ao and ¢2. 
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and then 
f(r) = exp (— Y) + ao exp (—X) in exp {—oxr — fonrr®}rdr 
0 


f(r) = exp (— Y) + exp (—X){1 — exp (-Y[4+ r"))} 
+ exp (—X + $Y/n’)(2Y/n)? 
x ferf (Y2L1 + m)) — erf (Y?/m)}. (41) 


Figure 15 shows curves of f(r). The ordinate Y? = (o7)'r was used 
as a convenient normalized range. It may be interpreted as the square 
root of the mean number of other peaks within range r of P;. As (40) 
shows, the parameter #(Z,) enters into f(r) through the variable X. 
The f(r) curves for different values of F(Z») lie close together at short 
ranges. As the range increases, f(r) falls more sharply for small E(2»). 
There is a limiting curve, as E(2») ~~, which is obtained from (41) 
by setting X = 0. As (40) shows, X = 0 also corresponds to the limit 
a—;i.e., this limit represents the planar model. 


EXACT (PLANAR) 
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Fig. 15—Probability f(r) that a random point at ground level is visible from a peak 
r miles away. 
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Table Vil— Mean coverage area C 


E(20) oC 
25 8.3 
50 12.2 

100 17.3 
200 23.5 
400 30.6 
800 38.6 
1600 47.4 


The approximations which led to (41) are poor when r is small. 
Fortunately, the planar model is good for smaller r. The curve labeled 
“exact (planar)” in Fig. 15 was obtained by a numerical integration, 
using an exact equation (89) for the planar model. The planar curve 
crosses the 0.5 probability level at Y? = 2.8. Then r = 1.307? is the 
range at which the odds of finding a clear path are even. 

The behavior of f(r) for large r can be obtained by replacing the 
error functions in (41) by their asymptotic expansions. The leading 
terms are 

f(r) ~ exp (— Y) + w*Y" exp (--X). (42) 


In Fig. 14, the f(r) curve starts to depart from the limiting curve at 
values of Y near H(2). For larger Y, the factor exp (—X) in (42) 
becomes small rapidly. In the planar model, exp (—X) = 1 for all Y; 
(42) then shows that f(r) ~ 7?/Y = r/(or?). 

To good approximation, the integral (36) for the mean coverage 
area can be replaced by 


Cz i: " f(r)d(ar?) =o i ” f(aY. (43) 


The main contribution to C in (36) comes from the range 0 Sr 
< (4m — v)a, in which (89) holds exactly and (41) approximately. 
Then (41) will be used for f(r) in (48) and, since f(r) — 0 rapidly for 
larger r, the range of integration has been extended to0 S$ Y < o~. 

Since (40) and (41) express f(r) in terms of Y and the single param- 
eter H'(2o), (48) shows that oC is a function of H(2o) only. Table VII 
gives values of this function, obtained from (41) and (43) by numerical 
integration. These values also represent mean numbers of peaks 
visible from a random point on the ground, as (37) showed. 

In Table VII, oC appears to be a slowly increasing function of 
E (20). As H(2o) —~, the model becomes planar and then f(r) ~ 7?/Y. 
Then (43) shows oC —~ in the planar model limit. 
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Two Design Techniques for Digital 
Phase Networks 


By F. J. BROPHY and A. C. SALAZAR 
(Manuscript received October 7, 1974) 


Two computer-aided algorithms for the design of all-pass digital filters 
are presented. The first technique is based on a linear programming ap- 
proach to solving the approximation problem posed by the minimax design 
of an all-pass digital filter. A new iterative algorithm with stability con- 
straints is offered for direct form design. The second technique implements 
a gradient search for those quadratic factors of an all-pass transfer function 
that lead to a locally optimal approximation (in the least-squares sense) 
of a desired phase function. New initial guess procedures and the parame- 
terization of linear-phase offset enhance the least-squares design procedure. 
Examples illustrating the result of both procedures are included. 


I]. INTRODUCTION 


The increasing availability of digital signal processors such as those 
described in Refs. 1 and 2 has generated much interest in the algo- 
rithmic design of digital filters. One particular class of recursive 
digital filters commonly referred to as all-pass digital networks has 
an important and interesting design problem associated with it. That 
is, the design objective for this type of filter involves the following 
transfer function 
by_,27* 

& 


H(z) = 


Mz uM= 


bpe7* 


k=0 


Because of the relationship between numerator and denominator 
polynomials, the number of degrees of freedom in filter design has 
been reduced to N from the usual 2N. Since the magnitude function 
of H(z) is precisely 1.0 on the unit circle, the design problem is 
focused directly on the phase variation of H(z). The importance of 
this design problem does not arise from an academic viewpoint. 
There are signal processing applications in which an influential 
factor in signal fidelity is the amount of phase distortion present in 
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Fig. 1—(a) Original pulse. (b) Phase-distorted pulse. 


the medium. The effects of phase distortion in communication systems 
are illustrated in Refs. 3 and 4. Apart from nonlinear phase equalization 
applications, all-pass networks can be used to provide a constant 
phase shift over a specified frequency band or bands. The Hilbert trans- 
former commonly found in bandpass modulation systems is just one 
example of this application. In constructing phased arrays in radar 
and seismic research, constant phase shifters are also found to be 
useful.*.§ Figure 1 illustrates how a constant phase offset can shape 
(or distort) the impulse response of a system where f(t) and f*(t) 
differ by a constant phase offset of r/2. A constant phase shift of any 
amount besides an odd multiple of 1/2 will produce a pulse with a 
single large lobe. Equalization of this type of distortion is again 
possible by all-pass networks. 

Previous work’ has addressed the envelope delay design problem. 
In many cases, this is sufficient but, as seen above, there are applica- 
tions where the phase function must be treated directly. 

Our design techniques are for all-pass structures where the design 
criteria stem from the phase function directly. The first technique, 
described in Section II, is a new method for designing all-pass networks 
using linear programming. This approach allows for fast (at least 
quadratic), always convergent design of phase networks. For the 
first time, stability can be treated directly in the design procedure. 
The second algorithm is based on a gradient search procedure on a 
least-squares criterion. The basic approach is analogous to those 
described in Refs. 7 to 9. The all-pass structure reduces the number 
of variables and simplifies the gradient calculations. In addition to 
developing the algorithm, we provide initial guess procedures and 
linear-phase offset parameters that enhance the algorithm. These 
initial guess procedures are new noniterative filter designs that can 
serve as excellent all-pass approximations in their own right. 


ll. A LINEAR PROGRAMMING APPROACH 


A need for fast, reliable design of all-pass digital filters has been 
shown in the previous section. Linear programming techniques have 
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been found to be useful in rational function approximations!" and 
have been applied to the magnitude-only design of digital filters.!2.4 
Here we show how the all-pass structure in digital filters can be trans- 
formed into a problem that can be handled by linear programming 
techniques also. As we shall see, the rational function differs from the 
magnitude-only case. Most importantly, this technique allows the 
question of stability to be handled directly in the design procedure. 
Other techniques that consider the phase or envelope delay variation 
of the digital filter (see Refs. 7 and 9 and Section III of this paper) 
deal with stability with a more heuristic approach. 

To develop the linear programming design method, we first recall 
that the all-pass transfer function is 


P(e) ws by + bn—127! + by_22* + +++ + Doz 











Be) = OE) by bet + Fe o) 
P(z) as 2 (bye% + bye" + pein + bo) : (1b) 
Q(z) (Gye™ Uys ee a 
Hence, the phase function of (1) on the unit circle is 
P(e) | ae = 
¢ (# os he 29LQ (27) J] j21=1. (2) 


From (2) we note that the phase variation of H(z7!) is equivalent 
(modulo a constant multiplier and an N sample delay term) to the 
phase of Q(z~"). Henceforth, we address the problem of synthesizing 
Q(z“). The phase variation of Q(z) on the unit circle is 


ca 2 b; sin 2arkf 
¢[Q (2) ]| je.—1 = tan-? —————_ 
d 6, cos 2rkf 
k=0 
tan o[Q (27) ]| a= Imag [Q(e-#*7) V/Real [Q(e-2*/)]. (3) 


Further, 
N 


> 0b, sin 2rkf 
tan ¢(Q(e-#"7)] = Sea s aC (4) 
2, o % COS 2akf 


Our design criterion is chosen to be 





R(fn) 
D(fn) — SM ty oy roy ’ 
igre (fn) S(fa) n=0,1,2 M 


where D(f) is the tangent desired phase function and M is a number 
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of frequency points* (>>) chosen to ensure adequate approximation 
over a subinterval of |f| < 4, namely, 0 < fo <fi:-- <fu < 34. We 
recall that the desired phase function has been scaled down by —4 
because of the factor of —2 appearing in (2) and will have a delay of N 
samples inherent in its design by the z— factor of (1b). It is important 
to note here that the norm is applied to the tangent of the desired 
phase function instead of the desired phase function itself. 

If we prevent S(f) from assuming the value zero, we seek the 
minimum value of A, 


Using the differential correction idea of Ref. 10, we expand the right- 
hand side of (5) in an iterative form: 


AS(fn) & ArSu(fn) + (A — Ax)Sk(fn) + BS (fn) — Su(fn) Jax. (6) 


The intention is to iterate toward those values of {b;} that minimize A. 
The subscript k indicates the kth iteration. We then have, from (5) 
and (6), 


|D(fn)S(fn) — R(fn)| — AwS (fn) — (A — Ar)Si(fn) S 0, 
which translates into a familiar pair of equations! 
[D(fn) + Ar IS(fn) — R(fn) + (A — An)Si(fn) 20 (7) 
[—D(fn) + ArIS(fn) + R(fn) + (A — A)Si(fn) 20. (8) 
Substituting the series forms for R(f,) and S(f»), we have 


N 


LX (LDF) + Ax] cos 2ajfn + sin 2ajfn}b; + (A — Ax) Se(Fn) 
= — D(fn) — Ar (9) 


zy {[—D(fx) + Ax] cos Qmjfn — sin Qmjfn}d; + (A — Az) Se(fa) 
a D(fn) pat Ai, (10) 


where bo = 1 is the normalization made. We have in (9) and (10) an 
over-determined set of 2M equations in N + 1 variables. The objective 
is to minimize A, one of the variables. It would seem that the condition 
S(fn) > 0 would be necessary to solve (9) and (10). But the phe- 


* An extension into a weighted criterion can be handled, but is suppressed in this 
presentation. M was chosen to be in the range 4N (N large) S$ M S 10N (N small) 
in our implementation of the algorithm. 

Therefore, the nonlinear nature of the tangent transformation may inhibit 
designs in the neighborhood of z. 
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nomenon experienced in Ref. 10 occurs here also. That is, if Si(fn) > 0, 
0snsM, then Siii(f,) > 0,0 Sn S M, also. 

Standard linear programming techniques can now be used on (9) 
and (10) to iterate toward a minimum A. However, no restriction has 
been made on the locations of the zeros of Q(z~!). Now there exist 
sufficient conditions for stability that can be written as linear con- 
straints. We have looked at two of these, e.g., a restriction that 
bi, be, ---, by of (1) form a monotonic sequence" or the restriction that 
the sum > 2-0 b; cos 2rkf > 0, wfE[0, 4].15 (The formulation of the 
linear programming problem gives us this condition on the subset 
of [0, 4] over which we are approximating.) For an example of a filter 
designed using this technique and the latter constraint to assure 
stability, refer to Fig. 2. Curve B is the sixth-order approximation to 
Curve A (only approximated over [0.075, 0.425 ]). 

However, the filter designer may decide that these types of con- 
straints are too restrictive for his particular applications. Nonlinear 
stability constraints, such as those found in Ref. 14, Chapter 3, can 
be included via the cutting planes algorithm,!* but this may require 
excessive computation times. Another suggestion involves interrupting 
the standard simplex method for solving the linear programming 
problem after each iteration. We may then further constrain the b 
vector used in the next basic feasible solution to a choice (i.e., some 
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Fig. 2—Sixth-order approximation using linear programming method. 
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Fig. 3—Phase error vs bandwidth for various orders of Hilbert transformers. 


“maximum’’) from among those vectors that would result in a stable 
filter in addition to the normal improvement of an object function. 
Using the standard formulation of the problem with no additional 
constraints or techniques necessary to assure stability, we were able 
to design many Hilbert transformer filters.* Figure 3 shows the relation- 
ship between the maximum error (recall that the tangent of the 
desired function is approximated) and a bandwidth (the filters were 


“rir designs of Hilbert transformers are well documented (see Ref. 17). There, 
90-degree phase is guaranteed, and the magnitude of 1.0 is approximated. 
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Fig. 4—Tenth-order Hilbert transformer. 


designed* over [f, 0.5 — f], f = 0.075, 0.05, ---, 0.225) for various 
orders of filters N = 4, 6, 8, 10. The log of the maximum error is given 
in the figure. 

The minimax approximation formulated here is performed on the 
tangent of the desired phase and not on the desired phase itself. For 
very good approximations, however, no penalties seem apparent. We 
have briefly looked at methods to design minimax phase approxima- 
tions based on the algorithm we have presented here. Our conclusions 
are that a two-stage design algorithm is required to iteratively locate 
a proper weight function that will “prewarp” the “tangent’’ design 
so that the weighted “tangent” design is minimax and the phase 
approximation is itself equiripple. 

Figure 4 illustrates the effect of the tangent transformation in this 
design procedure. In this figure, we see the phase of the resultant 
design (and its error function). This is a 10th order approximation to a 
90-degree phase shift over [0.05, 0.45]. While the design guarantees 
a minimax solution (equiripple) to the tangent, we can see that the 
resultant phase approximation is not exactly equiripple.t We have not 


* Each design only required a few (e.g., 5) iterations. 
We can see from Fig. 2 that the effect that the tangent transformation has on 
the error curve also depends on the values of the desired function. 
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implemented an algorithm to find the minimax solution to the phase, 
since, for our needs, the improvement in the phase approximation from 
the method outlined here did not seem to justify the use of a modified 
algorithm. 


lil. A GRADIENT SEARCH TECHNIQUE FOR LEAST-SQUARES DESIGN 


The next design algorithm we describe involves the computation of 
the gradient vector relative to the set of coefficients in a product of 
quadratic factors. The transfer function of an all-pass digital filter, 
expressed as a product of second-order sections, is: 


». ff btewtt+e2 
Hoe) = Uh (Pea): se 
The least-squares form 
L 
B= ¥ [D(fx) — Ang H(c-#*) Po(fi) (12) 


will be used as a measure of the approximation error from the desired 
function D(f) on the set of frequency samples {f,}{/. Here, w(f) 
denotes a nonnegative weighting function. A. G. Deczky has also 
considered gradient techniques applied to the least-square design of 
all-pass digital filters.’ In that paper, the emphasis was on envelope 
delay design. However, as shown in Section I, there are applications 
where envelope delay approximations are not adequate. Specifically, 
there are cases where phase distortion (e.g., constant phase offset) 
must be eliminated with an all-pass structure. 

Our design algorithm stems from familiarity with Ref. 8, which 
considers magnitude-only designs. With the least-squares criterion, 
the cascade second-order section form can be used. The advantage 
is that coefficient accuracy problems are minimized. As an alternative 
to the linear programming approach considered in Section II, this 
least-squares approach also enables one to more easily control the 
linear-phase offset permitted in the design. However, a disadvantage 
of the least-squares approach is that stability of the designed filter 
cannot be handled directly. Stability is obtained by confining the 
gradient movement to within the unit circle. This constraint may 
increase the likelihood of reaching an unsatisfactory local optimum. As 
we see later, there are initial guess procedures that provide excellent 
approximations to the desired phase function which, through the design 
algorithm, increase the likelihood of reaching a satisfactory local 
optimum. 
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3.1 Gradient calculations 


We find the entries of the gradient vector are 


oH do Ang aMe in fk) 


= —2 = [D( fx) — Ang H(e~?*4*) }w( fx) (13) 


dex; 
a =—2 > [D(f,) — Ang H (ee) w(fi) dang) (14) 


Here we define ¢(f) = Ang H(e~”"s) = tan“/(f)/R(f), where 
I(f) = Imag H(e-”"/) and R(f) = Real H(e-”*/). We seek 


29() = R(f)Ia(f) — I(f)Ralf) 
"ee = R(DIa(S) — LN Ra), 


where prime (’) denotes the partial derivative relative to the subscript. 
After some algebra, we find 


ot 2(1—B)FAf)snQsf isis M (15) 


5 = 2F(f)(sin 4af tasin2ef) 15i5M, (16) 


where F;(f) = |1 + ae~#"4 + B,e~#*7|-?. Finally, (13) and (14) can 
be simplified for 1 $7 S$ M to 


_ = — 41 — Bi) db CD (fi) — (fn) IF i(Se)w( fi) sin 2af; (17) 
oF 4 e [D(fe) — (fe) Fife) w (fe) (sin 42f, + a; sin 27f,). (18) 


The minimization of £ in (12) then proceeds with an iterative algorithm 
that is based on the formula 


cf) = cl) — 6, 1An( VE) a-1, (19) 


where c is the coefficient vector (a1, 81, a2, B2,-++,am, Bu) at the 
nth iteration, €, 18 the nth step size in the coefficient adjustment, A, 
is a positive definite matrix at the nth step (=/ in the case of the 
steepest descent algorithm) and (V£), is the gradient vector whose 
entries are given by (17), (18) at the nth iteration (we use the Fletcher- 


DIGITAL PHASE NETWORKS 775 


Powell algorithm). An initial guess procedure is required to start an 
iterative algorithm such as that of (19). 


3.2 Initial guess procedures for all-pass networks 


Convergence to a local minimum at which the approximation to a 
desired phase function is satisfactory can be made easier if a good 
initial guess is provided to c™ of (19). A desired feature of an initial 
guess procedure is that it be simple in nature. After all, excessive 
computation and effort should not be expected in simply starting a 
complex algorithm. In this section, we consider two procedures in 
which only a linear set of equations need be solved to obtain initial 
values for {b,}#=0 of (1). The value of having several initial guess 
procedures is that the designer may want to exercise his algorithm 
from multiple starting points to choose the best from a set of local 
optima. The following initial guess procedures operate on the direct 
form of H(z~)(1) which can be factored to the cascade form (11). 


3.2.1 Tangent approximation by Gauss’ method 


From (4) in Section II we know that a desired phase function can 
be approximated by considering a monotonic function of the phase, 
namely the tangent. Hence, 


b;, sin 2arkf 


i 


ol iM= 


tan ¢(f) = — (20) 


> 0b; cos 2xkf 
k=0 

is the approximating function of the tangent of half the desired phase. 
If we require the estimates of the desired phase tangent [tan ¢a(f) | to 
be “good” at a number of frequencies, we then have the following 
equations: 


N N 
tan da(fo) >> bg cos 2afyo — D> by sin 2rkfy 
k=0 k=1 


= 19 
N N 
tan da( fi) > b; COs Qrkfi _ a b; sin Qrkfi =] 
k=0 k=1 
N N 
tan ¢@a(fu) >> b, cos 2akfuy — > b, sin 2rkfy = ry. (21) 
k=0 k=1 


If {rn}° were all zero, then the approximation would be exact. The 
objective then is to minimize >>74) rz, where M > N. This problem 
is a least-squares minimization problem for which the solution is 
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derived from solving a set of normal equations: 


(a1, 41) (a1, a2)- ++ (a1, ay) | [Ox (a1, d) | 
(ae, a1) (a2, az)+ +> (a2, aw) be = (az, d) (22) 
(ay, ax) (ay, av)} (bv) (aw, 4) 
or 
Ab =e, 
where 


an = (tan ¢a(fo) cos 2rnfy — sin 2rnfo, 
-tan ¢a(f1) cos 2rnf1 — sin 2rnfi, ---, 
-tan da(fu) cos 2mnfy — sin 2rnfy) n= 1,2,---,N 


and 
d = [— tan ¢a(fo), — tan da(fi), ---, — tan ¢a(far) ] 
b= (bi, be, mae aes | by) and e = [ (a1, d), rane | (ay, d) J. 
Let 
= * Ss fe) 
Pmax — pene { Irs} and p= uM’ 


where r* = (79, 71, °°", Tu), the residual values after the least-squares 
approximation. If pmax — p is large (it is always positive), then a 
Chebyshev approximation may be desirable.'8 


3.2.2 Tangent approximation in Chebyshev sense 


It is well known that the minimax solution to (21) requires solving 
an appropriate subsystem of N + 1 equations. Further, the minimax 
solution of N + 1 inconsistent equations can be effected by examining 
the least-squares solution to the same set of equations and proceeding 
to solve a set of N linear equations. 

For our purposes here, an effective method of obtaining an initial 
guess for the iterative procedure implied by (19) is that of choosing 
M = N and obtaining the minimax solution to (21). This can be done 
by solving (22) for b = (b1, be, ---, by) and then evaluating (21) for 
the residuals 79, 71, -- +, ry. The minimax solution to (21) is then given 
by the linear set of equations 


Bb = a, (23) 
where B = (bj), N + 1 by N matrix with b;, = tan ¢a(f;) cos 2rkf; 
— sin 2rkf;, o = e[sign (ro), sign (71), ---, sign (rw) ], and 


cud *2 N * 
e= >> re/ > [rz|. 
k=0 k=0 
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It may be noted that only N of N + 1 equations are used in the solu- 
tion of (28). 


3.2.3 Discussion 


It should be noted that no constraint has been made on the initial 
guess procedures of A or B to ensure that the resulting digital filter 
is stable. In fact, if )>f-o b, cos k2af should ever change sign in |f| < 4 
or at least in the subinterval of approximation [fo, far], then the 
transition from (20) to (21) is not really valid since a division by zero 
is implied. Should }°?<o b; cos k2zf be strictly positive over |f| S 3, 
then stability results. (The interesting point is that stability can 
result even if the cosine series does change sign in |f| < 4). However, 
the point to remember is that the resulting initial guess may be 
unstable. In our experience, we have not encountered any serious 
problems using these initial guess procedures. 

We must further remark that the inherent N sample delay present 
in these approximations [see (2)] could present a problem when 
designing filters with M +N sample delays. However, we feel, 
intuitively, that since some delay is unavoidable, a delay of the order 
of the filter will not, for most applications, be overly restrictive. 

The last point to consider is that the initial guess procedure of Sec- 
tions 3.2.1 and 3.2.2 obtains a direct form estimate of the digital filter 
coefficients. What is really required for c of (19) are quadratic factors. 
We remark that we make the transition from the direct form estimate 
of (20) to quadratic factors by using a Bairstow quadratic factoriza- 
tion routine. 


3.3 Some considerations for least-squares design 


Often the engineering systems requirement of a digital filter can 
tolerate a linear-phase offset. While the systems engineer cannot 
always adapt to an arbitrary delay, there will usually be a range of 
delays permissible to him. How then can a designer incorporate these 
relaxations into the design mechanism? One technique for doing this 
is to add an acceptable delay to the desired function to create a new 
desired function and proceed from there. By designing filters for each 
of the permissible delays, one can choose from among the delays and 
their associated errors to decide which filter to implement. 

In Fig. 5 we can see the error function of a sixth-order filter* (B) 


* We have not tested the limit of the order of filters that can be designed by this 
method, but we have obtained a twentieth-order approximation (20-sample delay) 
to the desired function in this example. Quality initial guess procedures help us do 
this without excessive computation times. 
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Fig. 5—Error curves for initial and final sixth-order Hilbert transformer designs. 


designed for a delay of six samples. The desired function is a 90-degree 
phase shifting filter with the approximation having weight 1 on 
[0.08, 0.41] and 0 otherwise. Note the quality of the initial approxi- 
mation (A) using the first initial guess technique outlined in Section 
3.2. Of course, the disadvantage of presetting the delay is obvious; 
the choice of the optimal delay from those that are acceptable is not 
automatic but requires a separate design for each delay. However, 
eq. (12) can be expanded to include delay as parameter A 


He =D) (D(fz) — Ang H(e-¥**) — A2nfi Po(fr). 


An optimal A can be found analytically at each step in the gradient 
search and at convergence A will represent the amount of delay which, 
in conjunction with the filter, produces the best design.* Of course, we 
cannot expect that this delay will represent an integral number of 
samples or even a delay that the designer can tolerate. Figure 6 shows 
a desired function (A) (this curve is only shown where the weight of 
the approximation is nonzero) and its fourth-order approximation (B), 


*It is possible to include a constant phase angle as parameter B similar to the A 
used here. In such a case, our procedure becomes an envelope delay design technique. 
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Fig. 6—Fourth-order approximation to desired shape (A). 


which is by solving for optimal A. Allowing for arbitrary delay, the 
algorithm obtained this optimal design with a 4.1-sample delay. 

We can offer an heuristic solution to guarantee integer delays in an 
automatic fashion; namely, at each step (that is, at each calculation of 
A), the nearest acceptable delay* is used to replace A in the algorithm. 
This, of course, places a serious strain on optimality, although it does 
permit an automatic design procedure. 

As a footnote to this algorithm, we remark that there is a tendency, 
when working with procedures for designing filters in the cascade 
form, to use a previous optimal design of order 7 as the initial starting 
point in the design of filters of order n + 2. In the case of magnitude- 
only design, this is easily implemented since the appended second-order 
section can be initialized with magnitude 1. However, in the all-pass 
presentation there does not exist any second-order section that can be 
added which does not distort the overall phase when using a previous 
optimal design of order n to provide the initial guess for a design of 
order n + 2. And so the user of this algorithm must consider the effect 
of the appended second-order section if he does not want to obviate 
the value of a previous design toward providing an initial guess. 


* Nearest in the sense of greatest reduction of (12); “acceptable” here means 
‘Gnteger.” 
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The Effect of Small Phase Errors Upon 
Transmission Between Confocal 
Apertures 


By I. ANDERSON 


(Manuscript received November 14, 1974) 


The effect of small, periodic phase errors upon transmission between 
two coaxial, circularly symmetric apertures is considered when the aperture 
phase distributions are confocal and the amplitude distributions are 
gaussian. The results are applicable to loss calculations in beam wave- 
guide systems with imperfect lenses. When the periods of the phase errors 
are less than one-half the aperture radit, the total loss 1s approximately 
4(Bi + 63), where B1, B2 are the peak phase errors (in radians) on the 
apertures. Phase errors with periods greater than the aperture diameters 
are found to cause comparatively little transmission loss. 


l. INTRODUCTION 


The use of beam waveguide! systems for the transmission of informa- 
tion,” or for the transmission of power,’ necessitates the design of lenses 
(or cylindrical reflectors’) as focusing elements. In the design of these 
elements, it is desirable to estimate the degradation in performance 
caused by surface profile errors. Such degradation results in trans- 
mission loss and, in a communications system, will contribute to 
interference. Typically, the profile errors are associated with machining 
operations and, for lenses with circular symmetry, these errors are 
frequently circularly symmetric. The principal effect of the errors 
is to impart small, circularly symmetric phase perturbations to the 
field distribution adjacent to the lenses. The purpose of this paper 
is to calculate the reduction in transmission, caused by phase errors of 
this type, in a simple system comprising two coaxial, circular apertures 
as shown in Fig. 1. The field distributions on the apertures may repre- 
sent the fields in the aperture planes of two antennas or the fields on 
adjacent lenses in a beam waveguide system. 

In the absence of phase errors the transmission between coaxial 
apertures has been extensively studied by Kay,* Borgiotti,® Heurtley,’ 
and others, with the principal objective of determining that aperture 
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TRANSMITTING APERTURE RECEIVING APERTURE 
Ay A2 


Fig. 1—Coaxial, circular apertures with confocal field distributions &{(r;) and 
83 (72). 


field distribution which maximizes the transmission. When the aperture 
separation is greater than some minimum, it has been found that this 
field distribution corresponds to that of the lowest order mode in an 
open, confocal resonator. The phase distribution appropriate to this 
mode is obtained when the phase fronts on the apertures are confocal, 
i.e., are spherical with the center of curvature at the center of the 
other aperture. The appropriate amplitude distribution is well approxi- 
mated? by a gaussian curve. For this (optimum) distribution, the 
transmission between the apertures can attain surprisingly high values 
even when the apertures are separated by many aperture diameters 
(see Refs. 5, 6, and 7). Although the effect of periodic and random 
phase errors upon antenna gain and side lobe level has been investi- 
gated by several authors," little information appears to be available 
regarding transmission between two apertures when each has phase 
errors. In the case of transmission between two reflector antennas, 
Chu” has obtained an upper bound for the loss resulting from those 
phase errors produced by displacement of the feeds from the reflector 
foci. Yoneyama and Nishida!’ have considered transmission through 
a two-dimensional, confocal beam waveguide consisting of lenses with 
random phase errors. We compare their results to those of the present 
study later in this paper. 

In the following section the total transmission loss in a confocal 
system, with small phase errors on apertures with arbitrary ampli- 
tude distributions, is expressed in terms of the losses associated with 
each aperture when the other is free from phase errors. In the next 
section explicit expressions are derived, for two cases of practical 
interest, when the phase errors are sinusoidal and when the aperture 
amplitude distributions are gaussian. These expressions are then 
discussed and compared with results from the literature. 
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Il. TRANSMISSION BETWEEN CONFOCAL APERTURES WITH 

PHASE ERRORS 

Consider the circular apertures A1, Ae of radius a, a2 which are 
separated by a distance d > au, az, as in Fig. 1. It is assumed that the 
apertures are focused at each other such that the tangential electric 
fields in the apertures, when each is transmitting in the absence of the 
other, have the quadratic phase variation 


’2 


6i(r;) = E;(r;) exp (3 ee iy 4 = 1, 2, (1) 


where the r; are radial coordinates in the apertures. In the absence of 
phase errors the E;(r;) are real. If interaction is neglected, the trans- 
mission between the apertures is readily found’: from the results of 
Hu" and Kay®: 





1 1 1 2 
= D | i 7 Fyedridre ; (2) 
where 
Fy. => Ey (11) Ee(re)J o(nrire) rire (3) 
and 
1 1 x 
D=4 | fo ie@aienar [| Bars rad | (4) 


In these expressions the 7; are normalized so that r; = r;/a;. The 
Fresnel number n = kaia2/d, with k the wave number, and J» is the 
Bessel function of order zero. In the special case when the aperture 
separation is much greater than the aperture diameters, we see that 
n <1 and that the aperture phases are uniform, 1.e., Si(r;) = E;(r:). 
Substituting the small argument approximation for the Bessel function, 
x<1, Jo(z) 1, eq. (2) then reduces to the familiar Friis trans- 
mission formula}® 

AfA$ 








The effective aperture areas A{ are defined by 
f Ei(r)rdr| 
Ag = 2rai ——____—- 2 = 1, 2, (6) 


i. | 8i(r) [2rdr 
e.g., in the case of uniform illumination, 6{(r) = 1 and A{ = zai. 


The far-field transmission, 7 in (5), is also expressible in terms of the 
gains (G) of the apertures Ai and A2: 


Xr 2 
T= G0( 75), n<<& 1, (7) 
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where 
Ar A$ 


G= SH, 





i= 1,2. (8) 


Returning now to the discussion of phase errors, suppose that the 
phases in the apertures depart from the ideal (confocal) distributions 
by amounts ¢1(71) in Ai and ¢2(re) in Ae. The transmission 712 between 
the apertures in the presence of these phase errors is, from (2), 


1 1 fi ; 2 
Tr. = D | / J Fy. exp [j(¢1 + $2) ldridre} , (9) 
o Jo 


where the ¢,(r:) are abbreviated to ¢;. T12 is expressible as 
Ti2 = Ty) — AT x, (10) 


where 7'y is the transmission in the absence of phase errors and AT 12 
is the loss resulting from the phase errors. Let 


T; = Ty) — AT;, ~=1,2 (11) 
be the transmission between the apertures when there is a phase error 


on only the aperture A;. From Appendix A we then find 


AT. & ATi + AT. + R, (12) 
where 


1 1 1 1 1 
AT; = | f / Fyedridre i / Fy¢2dridre 
D 0 0 0 0 
1 1 2 
= | f / Pabadrdrs} ie z= 1, 2, (13) 
0 0 
2 1 1 1 1 
R= D | f J Fyedridre / / Fyodidedridre 
. 0 Jo 0 Jo 


1 fl 1 fl 
— i is Pipidridrs [ i Faxdsdradrs (14) 
0 


Equation (12) states that the total loss incurred by (small) phase 
errors ¢1(71) and ¢2(r2) on confocal apertures Ai, Ae is approximately 
equal to the sum of the losses associated with each aperture when the 
other is free of phase errors together with the term R. In the next 
section we evaluate the expressions for AT; and R when the aperture 
amplitude distributions are gaussian. 


Ill. GAUSSIAN APERTURES 


In the case of gaussian amplitude distributions, /;(r;) = exp (—air%). 
To simplify the analysis we assume the a; to be sufficiently large that 
the upper limits in the integrals may be extended to infinity. The 
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transmission 7’) in the absence of phase errors may now be obtained 
from (2), (38), and (4) by noting!* 


ic exp (—aer3)J o(nrire)redre = i exp (- vt) (15) 
2a 4as 


to give 
16n2a 19 


me (n? + 4aiae)?’ 


exp (—ai) <1, iS 12, (16) 
Apart from differences in notation, this expression is identical to the 
corresponding result obtained by Kogelnik!” for the coupling of 
(fundamental) gaussian modes. We find 


Tyo =1 when n = 2Varae, (17) 


ie., within the approximation exp (—a:) <1, there are optimum 
amplitude distributions that will ensure complete power transfer 
between the apertures for a given n. A detailed analysis,? or numerical 
integration, indicates this to be a satisfactory approximation when 
a; > 2.8, i= 1,2. For example, when ai = az = 2.36, the exact*® 
results for identical apertures are 7’) = 0.9931, n = 5.00, and the 
approximate results are 7’) = 1.00, n = 4.72. From (6), the effective 
aperture area, A‘, of a gaussian aperture is 
Pea ney (18) 
Qi 

As expected from physical considerations A{ decreases as a; increases, 
i.e., as the aperture field becomes more concentrated about the aperture 
center. 

In the case of transmission with circularly symmetric, periodic phase 
errors, we take 


(1%) amt Bi cos (yrs), i= 1, 2, 
where (19) 


Bi = ka, Y= 28s, 





8; is the peak value of the error in radians (with 6; the peak profile 
deviation) and 1; is the period of the error. For errors of period much 
greater than the circumference of the apertures, y; <1 and then, in 
(9), di(7:) & Bi, t = 1, 2. It follows, therefore, that Ty. = To, ie., 
to this order of approximation, small, slowly varying, circularly 
symmetric phase errors do not affect transmission between the aper- 
tures. In the general case of small errors, the transmission loss AT 12 
for gaussian apertures is found from (12) with (13) and (14) by 
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substituting for the #,(r,;) and ¢.(r;). From Appendix B, we have 


ADs _ peyifam(y'/2) — Dr) — VOD, 1-12, CO) 


r ; 4a; 
se Nn? + 4ajae (21) 


D(x) = exp (—2?) a exp (7?)dr (22) 





where 
and 


is the (tabulated) Dawson integral.!8 The index 7 = {7} when 7 = {3}. 
The term FR in (14) may be evaluated approximately in two cases 
of practical interest. In the first of these, the apertures are sufficiently 
far apart that n < 1 so that Jo(nrire) & 1 in (8). Fie is then separable 
in functions of r; and r2 and hence, from (14), R = 0. For this case, 


AT 12 & AT, + ATs, (23) 


i.e., the total loss is approximately the sum of the losses associated 
with each aperture when the other aperture is free of phase errors. 
The total loss is given by (20) and (23) in which y; simplifies to 


, Vi 
ix nel: 24 
Y oe ( ) 


It is noted from (7) and (11) that 


AT; _ AG; 
To Gy’ 








nK1, 1=1,2, (25) 


where G; = G; — AG; is the gain of aperture A; with the phase error 
¢:. Hence, (20) with (24) gives the fractional reduction in gain of the 
aperture A; resulting from a (small) periodic phase error. 

The second case of practical interest arises when the amplitude 
distributions on the apertures are optimized in accordance with (17) 
such that, in the absence of phase errors, the transmission is unity. 
From Appendix C, the term F in (12) is negligible in this case provided 
v1, ¥2>> nN = 2Varaz, i.e., the periods (J;) of the phase errors satisfy 


h< Ae and b<« xe ‘ (26) 
ae a1 


The transmission loss, AT12, between the apertures is then the sum of 
the losses associated with each aperture as given by (23). This result 


788 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1975 


implies that the transmission through a sequence of confocal lenses, 
each with small phase errors of period satisfying (26), may be obtained 
by calculating the transmissions associated with each lens in the absence 
of phase errors on the other lenses. Furthermore, (23) indicates that 
it is not possible to compensate for such phase errors on one lens by 
introducing phase variations on an adjacent lens. When (26) is satisfied, 
the Dawson integrals in (20) may be replaced by the first terms of the 
asymptotic expansion (44) to give 


AT; 483, i= 1,2. (27) 


It is of interest to note the physical significance of the condition 
(26) for the validity of the approximate forms (23) and (27). As 
expected from the theory of diffraction gratings, a circularly symmetric, 
periodic phase perturbation on a circular aperture generates!’ two 
additional side lobes in the aperture radiation pattern. If the period of 
the phase perturbation is J, then the two side lobes are symmetrically 
located about the main beam at an angle 6 = sin (d/l). Consider 
now the two apertures of Fig. 1 with phase errors of period 1; in 
A; and Il, in Ag. If the apertures are sufficiently far apart the side 
lobes, due to the phase error J; in Aj, will not intercept As provided 
sin— (A/l1) >> tan-! (a2/d), i.e., for small angles, li « Ad/az. Similarly, 
the main beam of A; will not couple energy to the side lobes of Ae 
provided 1. « \d/ai. The condition (26) implies, therefore, that energy 
is coupled from Ai to A» via the main beams alone. 

Figure 2 shows the transmission, as a function of y = y1 = y2, 
between two identical apertures as obtained by numerical integration 
of (9) with (19), and as obtained from the approximate result (23) 
with (27). The upper curve in the figure applies for a = 4, n = 8, 
8 = 0.36 and the lower curve for a = 2.36, n = 5, 6B = 0.18. In the 
absence of phase errors 7’) = 1 for these (optimum) distributions. The 
dashed lines correspond to the. approximation (23) with (27). As 
anticipated earlier, the transmission is essentially unaffected by phase 
errors of large period, e.g., when y S 7/2, 1.e.,1 2 2\d/a. The approxi- 
mate form (23) with (27) is seen to be within about 1 percent of the 
exact result when y = 2n, i.e., 1 S Ad/2a. As an illustrative example, 
consider a beam waveguide system of the type described by Arnaud 
and Ruscio? with \ = 3 X 10-%m, d 80m, a } 0.5m. The parameters 
of this system correspond to those of the lower curve in Fig. 2. Substi- 
tution shows that small, circularly symmetric phase errors of period 
1 = 2a on the lenses will cause negligible transmission loss, and that 
the approximation (23) with (27) is applicable for phase errors of 
period 1 S a/2. 
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Fig. 2—Dependence on y of phase error loss. 


IV. COMPARISON WITH PREVIOUS RESULTS 


To conclude this discussion on the effects of phase errors, we briefly 
compare the preceding results with the work of others. An expression 
for the gain (@’) of an aperture with small, periodic phase errors was 
given in (25). Consider the special case in which the period of the phase 
error is much less than the dimensions of the effective aperture, i.e., 
1< VA‘. From (18) and (19), this implies y >> Va so that the Dawson 
integrals in (20) with (24) may be replaced by the large argument form 
(44) to give, with (25), 


/ 

he (28) 
where G is the gain in the absence of phase errors. Since this result 
depends only upon the magnitude £ of the phase error, it is anticipated 
that it may apply to random phase errors with correlation lengths 
that are small compared with the dimensions of the effective aperture. 
Ruze! has examined the reduction in aperture gain caused by such 
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random errors and we compare (28) with his results. In the particular 
case of a sinusoidal surface error of rms value e on a parabolic reflector 
antenna, we have 8 = ki = 2v2ke. From (28), the gain with this small 


phase error is 
} 2 
@ wex | - (2)']. (29) 


This expression, derived here for a sinusoidal phase error, is identical 
to that obtained by Ruze in the case of a random error. As noted in 
Section I, Yoneyama and Nishida!* have examined the effect of random 
phase errors on lenses in a two-dimensional, confocal, beam waveguide 
system. Their approach is based on the concept of a statistical beam 
mode and this leads to a description of the field distribution, and 
transmission loss, in terms of an integral equation. A computer was 
used to solve the integral equation by numerical iteration from the 
solution in the absence of phase errors. It is interesting to find that 
the conclusions of their study, of a two-dimensional system with 
random errors, are similar to those obtained here for transmission 
between circular apertures with periodic phase errors. In particular, 
it was found that the transmission was not appreciably affected by 
phase errors with large correlation lengths and that the loss for a 
given error tended to a constant value for increasingly small corre- 
lation lengths. 


V. CONCLUSIONS 


We have examined the effect of small, periodic, radial phase errors 
upon transmission between two coaxial, circularly symmetric apertures 
with confocal phase distributions. Two cases of practical interest have 
been considered when the amplitude distributions on the apertures 
are gaussian. In the first of these the apertures are widely separated 
with phase errors of arbitrary period. The total loss is then the sum of 
the losses associated with each aperture and is given in terms of 
tabulated Dawson integrals. This result reduces to a known form 
when the periods of the phase errors are sufficiently small. The second 
case of interest applies to transmission through a beam waveguide 
system with imperfect lenses. When the periods (J;) of the phase errors 
on the apertures satisfy 1; S a:/2, (¢ = 1, 2), where a; is the aperture 
radius, the total loss resulting from phase errors is approximately 
2(6? + 6%), where 81, B2 are the peak phase errors in radians on the 
two apertures. A comparison, based on numerical integration, shows 
this to be within about 1 percent of the exact result in a typical case. 
Phase errors with periods 1; = 2a; (¢ = 1, 2) have comparatively little 
effect upon transmission. 
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APPENDIX A 
Derivation of (12) 


For small phase errors, (¢1 + ¢2) <1 and the exponential in (9) 
may be expanded to second order. Recalling that the H; are real, i.e., 
Fy, is real, we then find 


1 fl 2 
Tw x | f i Fyfl — $(@i + 6:)*andn| 


+ 1 f f Fy2(¢1 + 6s)dndral | (30) 


Expanding the first bracket and noting that 


1 1 fi 2 
D | f i Fi2(¢i + 6:)*dradra| S (61 + $2)4axT'o, (31) 


where (¢1 + ¢2)max IS the maximum value of (¢1 + ¢2), we have, to 


second order in ¢1, ¢2, 
Tr = Ty) — ATu, (32) 


where 


1 fl 1 fl 
AT 12 ee lf / Fydridre | / Fy2(¢1 + ¢2)?dridre 
o Jo o Jo 


= 1 f [ Fyo(gi + b:)drdrs}" | (33) 


Expanding the brackets gives (12) with (13) and (14). 
APPENDIX B 


Derivation of (20) 


Substituting £;(r:) = exp (—a:ri), o:(r:) = Bi cos (yrs) (¢ = 1, 2) 
into (13) with (8) and (4) gives, with (15), 


AT; 





To = 2nB3L11(y:) re 213 (yi) | t= 1, 2, (34) 
where 
Invi) = f° exp (—ar®) cos* (ysr)rar, (35) 
0 
I2(yi) = ye exp (— 7?) cos (y:r)rdr (36) 
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and 





1 
n= 7 (n? + 4arae), (37) 
with 7 = {?} when 7 = {3}. Expanding cos? (y:r) and integrating: 
1 
Ii(yi) = tn + $12(27:). (38) 
Integrating by parts, 
Tove) = ge [1 — re [exp (— a) sin (var |, (89) 
which is expressible!® in terms of the (tabulated) Dawson integral, i.e., 
ib 2928) 
hea == |i 40 
He) = 5, [1 Tew ( (40) 
where 
D(x) = exp (—2?) I exp (7°)dr. (41) 
0 


From (34), (38), and (40), we then obtain (20) and from (19) and (37) 
we obtain (21). 


APPENDIX C 
Approximate Evaluation of R in (14) 


We derive an approximate expression for R when yi, y2 > n. It is 
assumed that the amplitude distributions on the apertures are opti- 
mized such that 7) = 1 with n = 2Vo1a2, where a1, a = 2.3. Consider 
the integrals in (14): Extending the integration limits to infinity and 
substituting F(r;) = exp (—air), 7 = 1, 2 gives, with (8), (4), 
and (15), : 


1 1 1 
i, i Fyedridr2 = Ine ; D= 4n* (42) 


Similarly, substituting ¢; = 8; cos (y:r:;) and using (15) and (40), 


Ei _ Bi fy _ riV20. vm pet). 
i i Fyodidridre => On? E A D ( n 9 (43) 


Since 71 >> 7, the Dawson integral may be replaced by the first two 
terms of the asymptotic expansion,” 





1 © 13---(2m — 1) 
to give 
DP Rid andes 45) 
J, J, Putdndna re — a ( 
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Further, 


1 fl 1 
/ J Fyodigdedridre = B1B2 : exp (—airj)ri cos (yiti)S(71)dri, (46) 
o Jo 0 


where 


(71) = i . exp (—aers)J o(nrire)r2 COS (yar2)dre. (47) 
Substituting the integral representation of the Bessel function 

Jo(x) = if : cos (x sin 6)dé, (48) 

interchanging orders of integration and expanding the cosine product, 

s(n) = 5 [LI + 1(-0) Ja, (49) 


where 
I(@) = i exp (—aers)re cos [ (ye + nri sin 6)re |dre. (50) 
0 


From (40), 
ee E Seo) 
Jar Nets v2 1 


D | = (ye + nri sin ||. (51) 





Since ye > 7, both Dawson integrals in (49) may be replaced by the 
large argument form (44) to give 


S(71) S10) & — yr”. (52) 
Evaluating (46) by (40) and using (44), 





[ [ Fyedidedridre ~w a (53) 


Substituting for the integrals in (14) and reducing (20) then gives 
Rk ~ _ 8182 8n? 











AT: + AT,™~ = 62 + BB vid oe) 
But |8182/(6i + 83)| S 3, Le., 
2n \? 
— SG) ee 
Since y1y2 >> 2n, (55) is much less than unity and so, from (12), 
AT 12 % AT; + AT>. (56) 
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Toward a Group-Theoretic Proof of the 
Rearrangeability Theorem for 
Clos’ Network 


By V. E. BENES 
(Manuscript received October 8, 1974) 


Methods from group theory and combinatorics are used to prove the 
(Sleptan-Duguid) rearrangeability theorem for Clos’ three-stage network. 
The nr-permutations realizable in such a network can be represented as a 
product Go 1H 9G, where G, H are subgroups realized by stages and 
gy ts the special cross-connect field used in making frames. Thus, rearrange- 
ability can be cast as GgH ¢G = Sar = symmetric group of degree nr. 
Since it 1s an elementary theorem that a permutation group containing all 
transpositions is symmetric, it is enough to show that the product Ge H ¢G 
is closed under multiplication and contains all transpositions. We prove 
that closure of the product is equivalent to a property of suitable partitions: 
existence of systems of common representatives. This property, formulated 
by J. B. Kruskal, ts a consequence of Hall’s theorem on distinct representa- 
tives. It is easily seen that Gg H ¢G contains all transpositions, so the 
Slepian-Duguid theorem follows. 


1. INTRODUCTION 


In this paper we continue the exploration begun in previous work!-° 
of the relationships between permutation groups and connecting 
networks that are made of stages, frames, and cross-connect fields. 
Our results concern a well-known theoretical result of this area, the 
Slepian-Duguid theorem, which states that Clos’ three-stage network 
with square switches is rearrangeable, i.e., realizes any permutation. 
Since the permutations realizable by a stage form a special kind of 
subgroup, the theorem has been viewed in terms of group theory as a 
factorization of the symmetric group S,, of degree nr into a product of 
three subgroups or, alternatively, into a product of two mutually 
inverse double cosets.’ 

We further illuminate this basic rearrangeability theorem by giving 
it as nearly group-theoretic a proof as we have been able to find. This 
proof starts from the known characterization! of the nr-permutations 
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realizable by a Clos’ three-stage network as a product Gy H¢G, 
where G, H are subgroups realized by stages and ¢ is a “canonical” 
cross-connect field. It then shows that this product is closed under 
multiplication, and that it contains all nr-transpositions, whence 
immediately, by an elementary theorem, that it contains any nr- 
permutation, 1.e., that Snr = Ge AH ¢G. 

In the course of this proof we show that the basic combinatorial 
backbone of the rearrangeability theorem is really the existence of 
systems of common representatives (scrs) for pairs of partitions. 
Since, in apparent contrast, Duguid’s original proof used Hall’s 
theorem on systems of distinct representatives (sprs) of subsets, we 
have also sought to clarify just how the rearrangeability result depends 
on Hall’s theorem. The contrast above is apparent only because there 
are standard ways of proving scr results from spr results. In the 
present context, the two approaches are equivalent and lead to the 
same results. However, the scr formulation is closer to the group- 
theoretic aspects than is Duguid’s original spr proof: it provides an 
scr property that is a consequence of Hall’s theorem and is necessary 
and sufficient for the product Gy—H¢G to be closed. The property 
was first formulated by J. B. Kruskal in unpublished notes about 
rearrangeable networks dating from 1964. 


Il. SETTING AND FORMULATION 


We now sketch the group-theoretic interpretation of the Slepian- 
Duguid theorem in some detail, as has been done in earlier work.® 
Figure 1 shows Clos’ three-stage network, composed of three sym- 


CROSS—CONNECT 
FIELD y-! 





LAST STAGE MIDDLE STAGE FIRST STAGE 


Fig. 1—Ge"'H¢G describes the permutations realizable by Clos’ three-stage 
network. 
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Fig. 2—Getting transpositions in Gg 4H ¢G: terminals on different outer switches. 


metrically placed stages interconnected by the ‘“‘canonical”’ cross- 
connect field ¢ and its inverse. Each stage can realize precisely those 
permutations from a certain subgroup of S,,, depending on the size 
and number of switches in the stage. The r n X n switches of each 
outer stage realize a subgroup G isomorphic to (S,)’, viz., all those 
that permute the sets {kn + 1,kn + 2,---, (kK +1)n}, k=0,---, 
r — 1, within themselves. A similar statement holds for the center 
stage, but with n and r interchanged, to define a subgroup H isomorphic 
to (S,)*. . 

Thus, if we think of the network in Fig. 1 as acting from right to 
left, and if we interpret composition of permutations as left-multipli- 
cation of the inner permutation by the outer, then the permutations 
realizable by Clos’ three-stage network with square switches are 
- precisely those in the complex 


Go FH oG. 


The Slepian-Duguid theorem says that this complex is exactly the 
symmetric group Sn, of degree nr. We note for future reference that 
all transpositions are realizable; this can be seen from Figs. 2 and 38, 
in which the remaining terminals (not shown) are connected through 
to ‘themselves,’ as is possible and indeed necessary to realize a 
transposition. 


NN Ps 
~ got 
s s 
e ~ e P e 
SA wo 
e xX e@ os e 
~~ -* 
e ~ = e s e 
_ Pe 


Fig. 3—Getting transpositions in Gg" yG: terminals on same outer switch. 
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lil. SYSTEMS OF DISTINCT REPRESENTATIVES 


Let X be a set, and Xi, ---, X» finite subsets of X. We make the 
following definition. 


Definition 1: Elements 11, +++, &m from X form a system of distinct 
representatives (SDR) of X1,:-+, Xn tffxui C X: anda; ¥ 2; if 1 ¥ Jj, for 
i,j = 1, +++, ™. 


Hall’s theorem$ gives a necessary and sufficient condition for the X; 
to have an spr, thus: 


Theorem 1 (Hall): Xi, +++, Xm have an spr iff for k = 1, ---, m, the 
union of any k X; has at least k elements. 


This result was used by Duguid in his proof of the rearrangeability 
of Clos’ network with square switches. It enabled him to decompose 
any permutation into a union of submaps each of which, in switching 
terminology, carried exactly one terminal on each input switch onto 
images that were spread over all the output switches. These submaps 
could then be accommodated, one each on a middle switch. 


Iv. SYSTEMS OF COMMON REPRESENTATIVES 
Let P = {P;} andQ = {Q;} be partitions of aset X with |P| = |Q|. 
Definition 2: A subset HE C X is called a system of common representa- 
tives (scr) for P and Q iff 
(AO) Ps} = 1, P,EP 
Ryser® gives an spk argument to prove a necessary and sufficient 
condition for two partitions as above to have an scr. In the cases of 


interest to us here, a sufficient condition can be given in a particularly 
simple way. We make 


I 


Definition 3: Q is an (r, n)-partition iff |Q| =r, and |Q;| = n for 
Q: € Q. An (1, n)-partition of X ts one into r sets each having n elements. 


We use substantially Ryser’s argument® to prove the following 
special case (Theorem 2.2, p. 51, of Ref. 5) of his result: 


Theorem 2: Let P, Q be (r, n)-partitions of X. Then P and Q have an scr. 


Proof: For j = 1, ---,r, let A; = {7:P; meets Q;}. Take any union of 
k, of these sets, A;,U --- U A;,, and observe that Q;,U --- U Q;, 
has precisely nk elements in it. Hence, at most r — k integers in the 
range 1, ---, r fail to be in some A;,, ---, A;,. Thus, 


|A;,U --+ U Aj, | = k, 
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so, by Hall’s theorem, {A,;} has an spr {7;}, and P;; (1 Q; ~ ¢. Hence, 
P and Q have an scr. 


V. ORTHOGONAL PARTITIONS 


We now prove a property of partitions that will later turn out to be 
equivalent to the closure of the permutations realizable by Clos’ 
network. 


Definition 3: Partitions P, Q are orthogonal, written P 1 Q, iff P: EC P 
and Q; € Q imply |P: (1) Q;| = 1. 
Remark: If P 1 Q, and x is a permutation, then tP 1 7Q. 

The next result was first given by J. B. Kruskal. 


Theorem 8: If P, R are both (r, n)-partitions, then there is an (n, r)- 
partition Q orthogonal to each of P and R. 


Proof: By Theorem 2, P and R have an scr Q:. Remove all elements of 
Q: from the P; and the Q; to give new (r, n — 1)-partitions P’ and Q’. 
Repeat to find Qo, Qs, ---, Qn, and then take Q = {Q;}. 

It is convenient to have notations for three special partitions which 
arise naturally from the switching applications we are making. Clearly, 
the inlets (or outlets) of the network in Fig. 1 can be partitioned 
according to what last (or first) stage switch they are on. Similarly, 
the “‘wires’’ of the cross-connect fields between the stages can be 
partitioned according to what middle switch they impinge on. Accord- 
ingly, we define the (r, n)-partition S (by ‘‘outer’” switches) as 


S= {S37 = 1, area e S; = {k:(j —1)n <k S jn}, 
and the (n, r)-partition M (by “middle” switches) as 
M = {M;,j = 1, ---, n}, M; = (hij = 1)r <k SS 9r}: 


It is also convenient to partition by terminal position on outer switches, 
so we define the (n, r)-partition T by T = {T;,7 = 1, ---,} with 


T; = {k:k =In+j forsome 0518 r- lI}. 
The canonical cross-connect field is defined by 


ena | 
n 





ej>i+| | + CG = 1) moa n3 7 = 1,2, ---, nr. 


The following properties can be verified: ¢T = M, S 1 T. Intuitively, 
g takes the jth terminal on the 7th switch into the 7th terminal in the 
jth switch. 
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VI. CHARACTERIZATION OF REALIZABLE PERMUTATIONS 


The next theorem will give a necessary and sufficient condition on a 
permutation 7 to be realizable in Clos’ network, i.e., to belong to 
Ge eG. We start with a lemma. 


Lemma: Let P be any (r, n)-partition. If there ts an (n, r)-partition R 
such that 
Pes; 


then there exists an element g © G such that 
ggP 1 M. 


The practical import of this result is as follows: Consider a frame of 
r n X n switches followed by n r Xr switches, with the canonical 
cross-connect field 9 in between (Fig. 4); then, under the hypothesis 
there is a setting of the right-hand switches (i.e., the rn X n), which 
has the effect of connecting each set of P to some terminal on every 
switch of the left-hand stage of nr X 1, i.e., it images each P,; so as to 
reach every left switch (exactly once). 


Proof of lemma: Let R = {R,}. Each R; is simultaneously an spr of P 
and one for S. Thus, if-we connect the terminals of R; to the first 
left-hand stage switch, we will have used up one terminal from each 
P-set and also one from each switch on the right. This procedure can 
be repeated with Re, Rs, ---, R, to give the result. Evidently, this set of 
connections defines an element g € G such that each set of ggP is 
spread over the left-hand stage switches, i.e., such that ggP L M. 


Theorem 4:7 € Ge oG iff there is an (n, r)-partition R such that 
SLR L ws. 


IMAGING OF P ONTO LEFT—HAND SWITCHES 


rxer nxn 


«~~ 

~PIS AN 
(r,n )—PARTITION OF 
THESE INLETS 





yg g 


TS EACH SET OF ygP 1S SPREAD 
OVER LEFT STAGE SWITCHES 


Fig. 4—Import of the lemma. 
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EACH SET OF yg3-1S EACH SET OF yg, 

IS SPREAD OVER THE IS SPREAD OVER THE 

MIDDLE SWITCHES ~ _7 MIDDLE SWITCHES 
~ 





93 ga! 92 Y 4 


T=93 9-' 92 p 94 


Fig. 5—ggiS 1 M 1 og;'8. 


Proof: Let M be the partition of nr by middle switches, i.e., the (n, r)- 
partition consisting of the n sets 


re Lgre 22g te J = O05 ke id, 


and note that hM = M for h C H. Suppose now that « € Ge oG 
with 7 = g3¢7!g2¢91 and gi, gs € G, and gz € H. It can be seen from 
Fig. 5 that each set of ¢g3'S is spread over all the middle switches. 
Similarly, each set of ¢g1S is spread over the middle switches. Combina- 
torially, and without the help of pictures, these facts follow from 
eT = M, from gS = S for g € G, and from S 1 T, and they can be 
rendered as 
egS 1 M 
gg3 'S L M. 


It follows from the observation above that g.M@ = M, and thus, by 
the remark after Definition 3, 


gggS LM 1 vgs'S, 


whence 
cS a gsp (1M + S 


or 
S 1 gr'g gs'M L r8. 


For R, we take gy‘ ¢~gz'M, and the necessity is proved. 
For the sufficiency, we use the lemma, according to which the 
hypothesis implies that there is an element gi € G such that 


gg 3S L M. 
Thus, in Fig. 5, by setting up g1 in the right-hand stage, we can connect, 
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for each 7 = 1, ---, , the terminals of z—4S,, one each to a middle 
switch. It remains to define g2 for the middle stage by collecting those 
destined for S1, Se, ---, and gs for the left-hand stage by distributing 
within each of the sets S;, Se, --- in the left-hand stage. This is done 
precisely as follows: Define g2 by switching a terminal / to third stage 
switch 7 iff 
LE ggir8;. 

It follows that g“geygir—S; = S;. Then define gs by switching, within 
each final switch, ¢—ge¢gir—t to 7. Then t = gs¢ 92991 € Go A ¢G, 
as was to be proved. 


Vil. CLOSURE AND FACTORIZATION 


Theorem 5: Ge"H@G is closed under multiplication iff, for any two 
(r, n)-partitions P, Q, there ts an (n, r)-partition R such that P 1 R 1 Q. 


Proof: Let P, Q be given (r, n)-partitions. If Gg ¢G is closed, then 
it is a group that contains all transpositions, and so equals S,,. Hence, 
there exist permutations 7; and 72 such that 


mS = P, rs 'S = Q. 


Since G oH ¢G is closed, it is clear that 271 belongs to it. By Theorem 
3, or by inspection of Fig. 5, with + = mem, we see there is a partition 
N such that 

S LN L (mm)8; 
that is, 

1S al wiN an ws 1S. 


For the requisite partition R, take 7,N, and the necessity is proved. 

For the sufficiency, let m1, m2 €& Ge "HoG, and let P = 718, 
Q = 7;'S. Then, by the hypothesis, there is an (n, r)-partition R 
such that 

PLRLQ; 
that is, 
mS LR L wz'8 
SL wi, RR cafe (me71) 1S. 


Hence, by Theorem 4, r2r1 € Ge 4H ¢G, and we have proved that 
Go ¢G is closed. 


Theorem 6 (Slepian-Duguid) : 
Snr = Go H oG. 
Proof: Immediate from Theorems 3 and 5, since the right-hand side 


contains all transpositions and is closed. 
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Vill. FURTHER PROBLEMS AND COMMENTS 


Since H is a group, it follows that gH ¢ is also a group, one conju- 
gate to H, and that the Slepian-Duguid theorem can be cast as a 
decomposition 

Sa = U GrG 


rCe-lHe 


into disjoint double cosets, similar to the classical Frobenius’ de- 
composition. It is tempting to expect some sort of connection with 
Frobenius’ theorem here. One can speculate, in particular, that there is 
a proof of the Slepian-Duguid theorem from Frobenius’, obtained by 
specializing the requisite cosets to those of the form GG with z in the 
conjugate y—!1H y, and showing that only these need be considered. 

In conversation, Richard Stanley has indicated that, in another 
problem, also concerned with showing that a certain set of generated 
permutations was all of S,,, he had used the known result that a 
primitive group containing a transposition is a symmetric group. His 
remark stimulated our original approach to a “group-theoretic” proof 
of the rearrangeability theorem: one easily shows that, if Go-1H ¢G is 
a group, then it is a primitive group containing a transposition; the 
problem then became to show that it was closed, a property that 
turned out to be equivalent to Kruskal’s orthogonal partitions result 
(Theorem 3). Since closure was by comparison difficult to prove, and 
since it became clear that Ge'H ¢G contains all transpositions, the 
simpler proof presented here could be used, making the original side 
trip via primitive groups gratuitous. Stanley’s idea, however, is still 
a possible proof method for other networks that lead to less trans- 
parent groups of realizable permutations. 
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